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Abstract 

This  thesis  involves  the  modeling  of  self-sustained  oscillations  in  the  flow  past  a 
rectangular  cavity.  The  emphasis  is  on  developing  low-dimensional  models  that 
are  suitable  for  analysis  using  tools  from  dynamical  systems  and  control  theory. 

Two-dimensional  direct  numerical  simulations  are  performed,  and  indicate  the 
presence  of  a  ‘Svake  mode,”  which  has  been  observed  previously  in  experiments,  but 
w^hich  is  much  less  well  understood  than  the  “shear-layer  mode”  usually  observed. 
We  characterize  the  flow  in  both  shear-layer  mode  and  wake  mode,  and  provide 
a  criterion  for  predicting  the  onset  of  wake  mode,  as  a  function  of  the  various 
geometrical  and  flow-related  parameters. 

We  focus  on  the  modeling  of  shear-layer  mode,  and  employ  two  distinct  mod¬ 
eling  approaches:  first,  we  use  the  method  of  Proper  Orthogonal  Decomposition 
(POD)  and  Galerkin  projection  to  reduce  the  Navier-Stokes  equations  to  a  low¬ 
dimensional  system  of  ordinary  differential  equations  (ODEs).  We  extend  the 
method  to  compressible  flow's,  using  approximations  that  are  valid  for  cold  flow^s 
at  moderate  Mach  number.  In  a  compressible  flow,  both  the  kinematic  and  ther- 
mod\Tiamic  variables  contribute  to  the  total  energy,  and  an  inner  product  is  intro¬ 
duced  w'hich  respects  this,  and  allows  one  to  use  vector-valued  POD  modes  for  the 
Galerkin  projection.  We  obtain  models  in  the  form  of  ODEs  with  between  2  and  60 
states,  and  compare  models  based  on  scalar-valued  and  vector-valued  POD  modes. 
All  of  the  models  work  well  for  short  times  (a  few'  periods  of  oscillation),  but  the 
models  based  on  scalar- valued  modes  deviate  for  longer  times,  while  in  general 
the  models  based  on  vector- valued  modes  retain  qualitatively  correct  djmamical 
behavior. 

In  the  second  modeling  approach,  we  model  the  underlying  physical  mecha¬ 
nisms  separately  (shear-layer  amplification,  acoustic  scattering,  acoustic  propaga¬ 
tion).  and  obtain  linear  models  that  are  suitable  for  control  design  and  analysis. 
We  design  a  controller  which  stabilizes  the  model,  and  implement  a  similar  control 
law  on  an  experiment,  demonstrating  significant  reduction  in  the  amplitude  of  the 
oscillations,  but  revealing  some  limitations  of  feedback  control. 
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Chapter  1 
Introduction 


1.1  Motivation 

Self-sustained  oscillations  occur  in  a  wide  variety  of  applications,  from  the  flow 
past  a  sunroof  in  an  automobile,  to  landing  gear  wells  and  weapons  bays  in  air¬ 
craft,  to  musical  instruments  such  as  flutes  and  organ  pipes.  In  all  of  these  flows, 
acoustics  play  an  important  role:  sometimes  the  acoustic  radiation  is  desirable, 
as  in  a  musical  instrument,  but  in  most  engineering  applications,  sound  produc¬ 
tion  is  undesirable,  as  in  the  car  sunroof,  or  aircraft  applications,  where  pressure 
fluctuations  can  be  so  intense  that  they  can  lead  to  fatigue  and  structural  failure 
of  the  aircraft.  Furthermore,  acoustics  often  play  an  integral  role  in  the  physical 
mechanisms  underlying  self-sustained  oscillations.  This  is  the  case  for  the  cavity 
flows  we  study  here. 

The  basic  configuration  of  a  cavity  flow  is  shown  in  figure  1.1.  The  main 
features  of  this  flow  are  the  formation  of  a  shear  layer,  which  amplifies  flow  dis¬ 
turbances.  and  the  subsequent  scattering  of  these  disturbances  into  acoustic  waves 
at  the  downstream  corner.  These  acoustic  wnves  propagate  upstream,  and  excite 
further  disturbances  in  the  shear  layer,  creating  a  feedback  loop,  which  often  leads 
to  self-sustained  oscillations  at  discrete  resonant  frequencies. 

In  the  incompressible  limit,  a  similar  mechanism  exists,  but  instead  of  pro¬ 
ducing  acoustic  waves,  vortical  disturbances  moving  past  the  downstream  corner 
produce  an  irrotational  field  w^hich  excites  further  shear  layer  disturbances.  This 
irrotational  field  is  analagous  to  the  acoustic  field  generated  at  higher  Mach  num¬ 
bers,  only  the  propagation  speed  is  effectively  infinite. 


Figure  1.1:  Basic  configuration  of  a  cavity  flow\ 
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1  Introduction 


The  focus  of  this  thesis  is  to  develop  tractable  dynamical  models  of  these 
phenomena,  motivated  by  the  desire  to  use  active  control  to  reduce  or  eliminate 
the  self-sustained  oscillations,  for  instance  by  moving  a  flap  or  injecting  mass  into 
the  flow  periodically.  Modem  control  theory  provides  many  powerful  tools  both 
for  designing  feedback  control  system.s,  and  for  analyzing  limitations  of  control. 
To  use  these  tools  most  effectively,  one  first  needs  a  model  of  the  system,  tjpically 
in  the  form  of  a  system  of  ordinary  differential  equations  (ODEs).  The  motion 
of  fluids  is  described  very  accurately  by  the  Navier-Stokes  equations,  which  are  a 
system  of  nonhnear  partial  differential  equations  (PDEs)  which  are  not  solvable 
analytically  for  any  but  the  simplest  flows.  To  apply  existing  methods  of  control 
theory,  it  is  necessary  first  to  reduce  the  complexity  of  the  system,  preferably  to 
a  low-dimensional  system  of  ODEs.  The  goal  is  not  to  describe  every  detail  of 
the  flow,  but  rather  to  eliminate  unimportant  details,  and  obtain  the  simplest 
possible  mathematical  model,  while  retaining  just  enough  of  the  details  to  be  able 
to  accurately  describe  the  general  flow  features  one  cares  about,  such  as  the  size 
of  the  oscillations,  and  the  deflection  of  the  shear  layer. 


1.2  Historical  perspective 

Historically,  the  study  of  self-sustained  oscillations  began  wdth  studies  on  edge 
tones  [Brown,  1937;  Powell,  1953,  1961],  in  which  a  jet  emerges  from  an  orifice 
and  impinges  on  a  sharp  edge.^  The  mechanism  for  edge  tones  is  similar  to  the 
mechanism  for  cavity  oscillations  described  above;  disturbances  are  amplified  by  a 
shear  layer,  and  scattered  into  acoustic  waves  by  the  sharp  edge;  the  acoustic  waves 
excite  further  instabihties,  closing  the  feedback  loop.  The  edge  tone  is  the  canonical 
problem  which  describes  sound  production  in  flue-type  musical  instruments  such 

as  flutes  and  organ  pipes,  and  is  a  clos:  relative  of  the  cavity  flow  discussed  in  this 
thesis. 

1.2.1  Cavity  oscillations 

Cavity  flows  have  been  studied  as  early  as  the  1950s,  with  the  measurements  of 
Roshko  [1955]  and  the  thesis  of  Krishnamurty  [1956],  who  points  out  the  connection 
with  edge  tones.  The  first  description  of  the  acoustic  resonance  phenomenon  gov¬ 
erning  cavity  oscillations  (illustrated  in  figure  1.1)  is  usually  attributed  to  Rossiter 
[1964].  who  provided  a  semi-empirical  formula  to  predict  the  frequencies  of  oscilla¬ 
tion  (this  formula  is  discussed  later,  in  chapter  3).  This  description  of  the  acoustic 
feedback  mechanism  was  known  earlier  for  edge  tones  (e.g.,  Powell  [1953,  1961]), 
but  Rossiter  was  the  first  to  explicitly  apply  this  mechanism  to  the  cavity  flow! 
and  provide  a  formula  for  the  resonant  frequencies.  Many  enhancements  to  the 
basic  theory  have  been  given  over  the  years,  and  are  described  in  chapter  3. 

Another  mode  of  cavity  oscillation  has  been  observed,  but  has  received  much 
less  attention,  and  is  relatively  poorly  understood.  In  incompressible  experiments 

(1961)  points  out,  an  early  description  of  self-sustained  oscillations  is  also  given  bv 
Lord  Rayleigh,  in  his  description  of  the  bird-call  [Rayleigh,  1945,  §371). 


1.2  Historical  perspective 
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for  an  axisymmetric  cavity,  Gharib  &  Roshko  [1987]  observed  a  wake  mode,  where 
the  oscillating  flow  over  the  cavity  resembles  the  wake  behind  a  bluff  body,  rather 
than  a  free  shear  layer.  Flow  features  in  this  wake  mode  were  qualitatively  very 
different  from  those  in  the  shear-layer  mode  described  by  Rossiter,  and  wake  mode 
was  accompanied  by  a  large  increase  in  drag.  Similar  dramatic  increases  in  drag 
had  been  previously  observed  by  Fox  [1968]  as  the  cavity  length  was  increased,  in 
flows  with  thin,  laminar  upstream  boundary  layers,  and  Roshko  [1955]  observed 
an  intermittency  “analagous  to  the  large  fluctuations  of  drag  which  occur  on  a 
bluff  cylinder  in  the  critical  range  of  Reynolds  number,”  where  the  flow  may  be 
switching  betv^een  shear-layer  mode  and  a  type  of  wake  mode.  Recent  experiments 
and  numerical  studies  by  Kriesels  et  al.  [1995]  for  a  flow  past  closed  branches  of  a 
pipe  (with  both  circular  and  2D  cross  sections)  also  demonstrate  flow  structures 
that  closely  resemble  wake  mode. 

More  details  about  cavity  oscillations  in  general  may  be  found  in  review  articles 
[e.g.,  Rockwell  &  Naudascher,  1978;  Blake  &  Powell,  1986;  Colonius,  2001];  see  also 
the  recent  analytical  work  by  Howe  [1997]  for  very  low  Mach  number  cavity  flows, 
and  Crighton  [1992]  for  edge  tones. 

Research  in  cavity  flows  has  been  renewed  in  recent  years,  largely  because  of  the 
possibility  of  using  active  control  to  reduce  the  amplitude  of  oscillations.  The  idea 
of  controlling  cavity  oscillations  is  not  new,  and  goes  back  to  Sarohia  &  Massier 
[1977]  and  Gharib  [1987].  Most  of  the  recent  approaches  to  control  have  been  either 
open-loop  (no  feedback)  [Vakili  Gauther,  1991;  Sarno  Franke,  1994;  Stanek 
et  al.,  2000],  or  have  involved  phase-locked  loops  or  adaptive  controllers  [Cattafesta 
et  al.,  1999;  Shaw  &  Northcraft,  1999;  Williams  &:  Fabris,  2000b]  which  do  not 
incorporate  models  of  the  cavity  physics.  The  results  of  active  control  have  not 
been  completely  successful:  often,  greater  benefits  are  obtained  by  using  passive 
techniques  such  as  fixed  spoilers  [Vakili  Gauther,  1991],  and  often  when  active 
control  is  used  to  reduce  the  amplitude  of  one  unstable  frequency,  the  amplitude 
increases  at  another  frequency  [Williams  &  Fabris,  2000b],  In  order  to  understand 
these  effects,  and  the  limitations  and  potential  benefits  of  active  control,  one  first 
needs  a  tractable  model  of  the  system.  Modeling  the  cavity  flow  is  the  focus  of 
this  thesis. 

1.2.2  Model  reduction 

In  order  to  obtain  low-order  models  for  the  cavity  flow,  we  use  some  tools  of  dynam¬ 
ical  systems  theory,  most  notably  the  method  of  Proper  Orthogonal  Decomposition 
(POD)  £Lnd  Galerkin  projection.  Galerkin  projection  is  a  method  for  obtaining  ap¬ 
proximations  to  a  high-dimensional  (even  infinite-dimensional)  dynamical  system 
by  projecting  the  dynamics  onto  a  low-dimensional  subspace.  Proper  orthogonal 
decomposition  is  a  method  which  uses  data  from  simulations  or  experiments  to 
determine  a  subspace  which  is  optimal  in  a  certain  sense. 

Also  known  as  Karhunen-Loeve  expansion,  the  POD  was  originally  developed 
in  the  context  of  probability  theory  [Loeve,  1978],  but  the  method  has  been  used 
subsequently  as  a  model-reduction  tool,  particularly  for  fluids  and  other  systems 
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described  by  PDEs.  In  the  fluids  context,  the  major  contributions  are  from  Lumley 
[1970],  Sirovich  [1987],  and  Aubry  et  al.  [1988],  who  used  the  method  primarily  to 
study  the  tiurbulent  boundary  layer.  More  recently,  the  method  has  been  used  to 
study  other  flows  (e.g.,  Khibnik  et  al.  [1998]),  but  so  far  all  of  the  flows  studied 
have  been  incompressible.  (Several  researchers  have  computed  POD  modes  for 
compressible  flows  [Ukeiley  et  al.,  2000],  but  have  used  the  POD  modes  only  to 
identify  coherent  structures  in  data,  not  to  develop  dynamical  models.)  One  of 
the  major  contributions  of  this  thesis  is  to  establish  a  framework  for  applying  the 
POD/Galerkin  method  to  compressible  flows. 


1.3  Overview  of  contributions 

The  main  contributions  of  this  thesis  are  as  follows: 

•  Direct  Numerical  Simulations  (DNS)  of  the  two-dimensional,  compressible, 
subsonic  flow  past  a  rectangular  cavity,  along  with  the  radiated  acoustic  field. 

•  Investigation  of  the  wake  mode  of  cavity  oscillations,  and  criteria  for  predict- 
ing  the  onset  of  wake  mode. 

•  Extension  of  the  POD/Galerkin  method  to  compressible  flows,  valid  for  cold 
flows  at  moderate  Mach  number. 

•  Comparison  of  reduced-order  models  obtained  from  scalar- valued  POD  modes 
and  vector-valued  POD  modes. 

•  Development  of  a  controller  to  reduce  cavity  oscillations,  based  on  linear 
conceptual  models  of  the  cavity  physics. 

•  Implementation  of  controllers  on  an  experiment,  demonstrating  noise  reduc¬ 
tion.  and  some  limitations  of  control. 

In  chapter  2.  wc  give  an  overview  of  linear  stability  analysis  of  parallel  flows, 
which  will  be  used  throughout  this  thesis.  Chapter  3  gives  a  description  of  the 
numerical  method  used  to  compute  the  flow  past  a  cavity,  and  describes  charac¬ 
teristics  of  the  flow  in  detail,  including  both  shear-layer  mode  and  wake  mode.  We 
propose  that  wake  mode  arises  when  the  amplitude  of  oscillations  exceeds  a  certain 
threshold,  possibly  creating  a  region  of  absolute  instability  in  the  shear  layer. 

Chapter  4  contains  the  model-reduction  methods  developed  in  this  thesis:  an 
overview  of  the  standard  POD/Galerkin  procedure  is  presented,  along  with  a  brief 
description  of  the  existing  methods  for  incompressible  flow.  One  of  the  main 
contributions  of  this  thesis  is  the  extension  of  these  methods  to  compressible  flows, 
presented  at  the  end  of  this  chapter. 

In  chapter  5,  we  apply  these  methods  to  the  cavity  flow,  and  compare  models 
obtained  from  scalar-valued  POD  modes  to  those  obtained  from  vector-valued 
modes. 


1.3  Oveniew  of  contributions 
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Chapter  6  presents  a  conceptual  modeling  approach,  based  on  the  Rossiter 
mechanism,  where  we  model  each  component  of  the  flow  separately  (e.g.,  shear- 
layer  amplification,  acoustic  propagation).  These  models  are  linear,  and  we  use 
them  to  design  controllers  using  the  Linear  Quadratic  Gaussian  (LQG)  design 
technique.  We  implement  these  controllers  on  an  experiment,  and  discuss  the 
results  in  chapter  6. 
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Chapter  2 

Linear  Stability  of  Parallel  Shear  Flows 


Since  the  late  19th  century,  linear  stability  theory  has  played  an  important  role  in 
modeling  of  shear  flows.  The  first  theoretical  contributions  appear  in  the  famous 
paper  of  von  Helmholtz  [1868].  who  established  the  appropriate  boundEiry  condi¬ 
tions  at  a  discontinuity  between  two  streams  with  dilfferent  velocities.  The  stability 
of  this  boundary  was  first  addressed  by  Kelvin  [1871],  who  looked  at  the  motion 
of  a  free  surface  of  water  as  wind  blows  over  it,  and  to  this  day  the  instability  of 
a  flow  with  two  parallel  streams  bears  the  name  Kelvin-Helmholtz  instability. 

The  foundations  of  the  modern  Unear  stability  theory  were  laid  by  Rayleigh 
[1880:,  who  used  his  new  framework  to  prove  remarkable  results  about  the  stabilit}' 
of  jet  profiles.  Rayleigh’s  theory  persists  today,  virtually  unchanged,  and  this 
theory,  with  some  modern  accoutrements,  forms  the  subject  of  the  present  chapter. 

Linear  stability  theory  will  be  used  frequently  in  subsequent  chapters,  to  de¬ 
scribe  the  behavior  of  the  shear  layer  spanning  a  rectangular  cavity.  The  de¬ 
velopments  here  wdll  be  from  this  perspective,  considering  the  appropriate  flow 
assumptions  and  boundary  conditions.  Mc^t  of  the  results  are  well  known,  but 
scattered  throughout  the  literature,  so  here  we  collect  them  in  one  place,  along 
with  a  few  recent  results,  such  as  the  adjoint  equations  and  biorthogonality  rela¬ 
tion  developed  in  sections  2.6  and  2.7,  We  describe  a  locally  parallel  theory,  for 
inviscid  flows.  We  develop  the  compressible  Rayleigh  equation,  along  with  the 
corresponding  adjoint  equations,  and  discuss  a  biorthogonality  condition  which  is 
useful  for  expanding  arbitrary  flow  disturbances  in  terms  of  eigenfunctions  of  the 
Rayleigh  problem. 

2.1  Equations  of  motion 

We  assume  that  the  flow  is  two-dimensional,  compressible,  and  isentropic,  and 
is  thus  described  by  the  Euler  equations  for  an  ideal  fluid,  written  in  Cartesian 
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coordinates  as 

(2.1) 
(2.2) 
(2.3) 

where  D/Dt  :=  5/(9t  +  u  - V  denotes  the  material  derivative.  The  density  p  and  ve¬ 
locity  u  =  {u,v)  are  nondimensionalized  by  the  freestream  quantities  poc  and  Uoo, 
the  pressure  p  is  nondimensionalized  by  PooU^,,  and  lengths  are  nondimension¬ 
alized  by  a  length  D,  as  yet  unspecified.  In  writing  the  energy'  equation  in  the 
form  (2.3),  we  assume  the  fiuid  is  an  ideal  gas  with  constant  specific  heats,  that 
there  is  no  heat  addition,  and  that  no  body  forces  act  on  the  fluid. 

Next,  we  Unearize  equations  (2.1)-(2.3)  about  a  parallel  mean  flow,  writing 


u{x,y,t)  =  U{y)  +  Su{x,y,t)  (2.4) 

y,  t)  =  6v{x,  y,  t)  (parallel  flow)  (2.5) 

P{x,y,t)  =  p{y)  +  Sp{x,y,t)  (2.6) 

p{x,y,t)  =  p+ 6p{x,y.t).  (2.7) 


We  obtain  the  linearized  equations  by  inserting  (2.4)-(2.7)  into  equations  (2.1)- 
(2.3)  and  neglecting  products  of  fluctuations: 


Dp 

Dt 


-f  p  div  u  =  0 


Dp 


Du  _ 

-f  7pdivu  =  0, 


Sp  -f  pSv  +  p  ( ~^Su  4-  —  )  =  0 

\9x  dy  ' 

Q 

Su  +  pU'Sv  +  —6p  =  0 
dx 

r  9  , 

6v  -I-  —6p  =  0 
dy  ^ 

,  \  f  d  ,  d  \ 


where  primes  denote  differentiation,  and  M  =  t/oc/aoc  is  the  free-stream  Mach 
number.  Here  we  have  used  p  ^  which  follows  from  the  equation  of 

state  for  an  ideal  gas  (a^  =  7Pdc/Pdc)-  Note  that  the  density  equation  decouples, 
so  we  need  only  consider  the  last  three  equations. 

Next,  w’e  write  the  fluctuations  as  normal  modes 


where  /  =  u,  i;,  or  p,  (2.8) 

where  q.o;  e  C.  Since  the  equations  are  linear,  a  linear  combination  of  solu¬ 
tions  is  also  a  solution,  so  by  superposing  normal  modes  we  may  construct  more 
complicated  solutions. 

B}  inserting  the  normal  modes  into  the  linearized  equations,  we  obtain  the 


2,2  Mean  profiles  and  boundary^  conditions 


9 


system  of  ordinary  differential  equations  (ODEs) 


where  ip  ==  and 


/  p{iuj  +  iaU)  pU' 


LuJ.Ol  — 


la 


p{iuj  +  iaU) 
A 

dy 


lOc 

d 

o  ^  - 


(2.9) 


(2.10) 


From  these  equations,  straightforward  Gaussian  elimination  yields  a  single  ODE 
in  p,  known  as  the  compressible  Rayleigh  equation: 


d^p 

dy^ 


2aU'  P'\^P 
u  +  at)  P  J  dy 


(a^  —  pM^{u  4-  ckU)"^)  p  =  0. 


(2.11) 


The  incompressible  Rayleigh  equation  is  recovered  by  setting  M  =  0,  p  =  1,  and 
is  usually  written  in  terms  of  t),  as  originally  derived  by  Rayleigh  [1880]: 


-  U%  =  0. 


(2.12) 


Note  that  both  of  these  forms  of  Rayleigh’s  equation  are  singular  when  U  =  —ujja. 
This  case  can  be  important,  and  is  addressed  in  section  2.5. 


2.2  Mean  profiles  and  boundary  conditions 

Later,  w^e  shall  use  the  above  equations  to  describe  the  shear  layer  over  a  rectangu¬ 
lar  cavity.  Here  we  discuss  the  mean  profiles  and  boundary  conditions  appropriate 
for  this  problem. 

For  the  shear  layer  spanning  the  cavity,  we  shall  use  two  different  velocity 
profiles:  actual  mean  profiles  measured  from  simulations,  and  a  hyperbolic  tangent 
profile  given  by 


e(j)  =  ^l+tanh|) 


(2.13) 


which  describes  a  free  shear  layer  with  (nondimensional)  speed  1  in  the  upper 
stream,  zero  velocity  in  the  lower  stream,  and  with  vorticity  thickness  5^,  Recall 
that  the  vorticity  thickness  of  a  shear  layer  is  defined  by 


U2-U1 

maxU'{y) 

y 


(2.14) 


where  U{y)  is  the  mean  velocity  profile  of  the  shear  layer,  and  U2  and  Ui  are 
the  speeds  in  the  upper  and  lower  streams,  respectively  [Brown  &  Roshko,  1974]. 
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(Here,  U2  =  1,  Ui  =  0.)  Another  quantity  frequently  used  is  the  momentum 
thickness 


0  = 


U{l-U)dy. 


(2.15) 


For  the  velocity  profile  (2.13),  the  two  are  related  by  6^  =  40.  As  we  discuss  in 
chapter  3,  the  vorticity  thickness  is  more  appropriate  for  describing  shear  layers 
over  cavities. 

Throughout,  we  will  assume  the  mean  density  is  uniform,  p  =  1.  This  assump¬ 
tion  is  good  for  cold  flows  at  the  moderate  subsonic  speeds  we  consider. 

The  developments  in  the  preceding  section  assume  that  the  flow  is  uniform  in 
the  streamwise  direction,  so  in  modeling  a  shear  layer  over  a  cavity,  we  necessarily 
neglect  the  presence  of  the  upstream  and  downstream  cavity  walls.  We  may  include 
the  effect  of  the  cavity  floor,  however,  by  imposing  a  wall  boundary  condition. 
We  denote  the  cavity  depth  by  D.  and  use  this  D  as  the  length  scale  in  the 
nondimensionalization.  The  wall  boundary  condition  is  then 


=  i.e.  ^  =  0^  aty  =  -l.  (2.16) 

The  remaining  boundary  condition  is  that  the  pressure  be  bounded  as  y  — ► 
oc.  Far  from  the  cav'ity,  the  mean  profiles  are  uniform,  and  so  the  governing 
equation  (2.11)  reduces  to 


^  -  {a^  -  +  af)p  =  0  (2.17) 

which  has  sol nt ion 

p  —  -r  (2.18) 

where  A“  =  -f  a)‘.  Since  the  pressure  must  remain  bounded  as  y  oc. 

this  implies  C2  =  0  (where  we  choose  the  standard  branch  of  the  square  root  for  A). 
The  farfield  boundary  condition  is  then 


p  =  C]€ 


as  y  ^  oc. 


(2.19) 


2.3  Method  of  solution:  spatial  vs.  temporal 

The  governing  equation  (2.11)  together  with  the  boundary  conditions  (2.16)  and 
(2.19)  determine  a  linear  homogeneous  boundary  value  problem,  and  hence  an 
eigenvalue  problem.*  Nontrivial  solutions  for  p  exist  only  for  certain  values  of  the 

*  Strictly  speaking,  in  an  eigenvalue  problem  the  parameter  must  appear  linearly,  w'hich  is  not 
true  for  either  u;  or  q  in  (2.11).  However,  both  u;  and  q  do  appear  linearly  in  (2.9),  so  this  equation 
should  be  thought  of  as  the  true  eigenvalue  problem,  and  (2.11)  a  convenient  representation  of  it. 


2.3  Method  of  solution:  spatial  vs,  temporal 
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parameters  u  and  a,  satisfying  a  particular  relation  of  the  form 


a)  =  0 


(2.20) 


for  some  function  F  (which  is  not  known  exphcitly).  This  equation  is  called  the  dis¬ 
persion  relation,  and  relates  frequencies  of  waves  to  their  corresponding  wavenum¬ 
bers. 

In  fluid  mechanics,  we  often  look  for  solutions  of  (2.20)  such  that  either  Im(a;)  = 
0  or  Im(a)  =  0.  Consider  the  former  case  first.  If  the  dispersion  relation  (2.20)  has 
a  solution  with  Im(u;)  =  0  and  Im(a)  <  0,  the  corresponding  normal  mode  (2.8) 
will  be  oscillating  at  a  particular  frequency  in  time,  but  exponentially  growing 
in  space.  Such  a  situation  is  called  spatial  instability.  Next  consider  solutions 
of  (2.20)  with  Im(a)  =  0.  If  there  is  a  solution  with  Im(u;)  <  0,  then  the  corre¬ 
sponding  normal  mode  will  be  sinusoidal  in  space,  and  exponentially  growing  in 
time.  This  situation  is  called  temporal  instability.  If  p  is  an  eigenfunction  with 
either  Im(a)  <  0  (with  E  R)  or  Im(a;)  <  0  (with  a  E  R),  we  call  p  an  un¬ 
stable  mode.  Correspondingly,  if  p  is  an  eigenfunction  with  either  Im(Q)  >  0  or 
Im(a;)  >  0,  we  call  p  a  stable  mode  (or  damped  mode).  If  both  CL;,a  E  R,  then  p  is 
a  neutrally  stable  mode.  We  consider  only  the  cases  when  at  least  one  of  uj,  a  is 
real. 

In  modeling  spatially  developing  shear  layers,  the  spatial  approach  is  more 
appropriate,  and  results  from  spatial  stability  analysis  agree  well  with  experi¬ 
ment  [Michalke,  1984],  so  it  is  the  approach  we  take  here.  The  general  procedure 
is.  given  a  frequency  cj  E  R,  to  look  for  a  value  of  a  E  C  such  that  the  boundary 
value  problem  (BVP)  given  by  (2.11),  (2.16),  (2.19)  has  a  solution.  This  procedure 
is  made  systematic  as  follows. 

We  start  by  choosing  an  initial  guess  for  q,  and  attempting  to  solve  the  linear 
two-point  BVP  by  shooting.  We  write  equation  (2.11)  in  first-order  form 


dp 

dy 


=  p 


pAf{uj  +  aUf)p  + 


2aU' 
u  +  aU 


T> 


with  boundary  conditions 

p(-l)  =  c  P(y  —  oc)  =  de~^^ 

p'(— 1)  =  0  p){y  — ►  oo)  = 


(2.21) 

(2.22) 


(2.23) 

(2.24) 


where  c  and  d  are  constants,  now  arbitrary.  To  solve  the  BVP,  we  integrate 
the  equation  from  y  =  “1  to  an  arbitrary  point  y  =  yo,  and  call  this  half  of  a 
solution  pi(y).  We  then  integrate  the  equation  from  y  =  oo  (in  a  real  computation, 
from  some  large  value  where  the  mean  fields  are  approximately  constant)  to  y  = 
yo,  and  call  this  part  of  the  solution  P2(y)-  For  the  BVP  to  have  a  continuous 
solution,  we  must  choose  the  constants  c  and  d  such  that  Pi(yo)  =  P2(yo)  and 
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p'liyo)  =  P2{yo)-  These  two  conditions  are  simultaneously  possible  only  when 
Pi/ Pi  =  P2/P2  at  yo,  or  when 

G{a)  =  {pip'2  -  P2p'i)\y^  =  0.  (2.25) 

Note  that,  holding  c,  d,  and  w  constant,  the  value  of  G  depends  only  on  a.  If 
G  is  anahtic,  then  we  may  use  Newtoh-Raphson  iteration  to  solve  the  algebraic 
equation  G{a)  =  0  for  the  eigenvalue  a.  The  normalization  of  the  eigenfunctions 
is  arbitrary,  and  we  choose  the  normalization  maxy  v  =  1. 

2.4  Number  of  unstable  eigenvalues 

For  typical  shear  profiles  U  that  we  use  in  this  thesis,  there  is  only  one  unsta¬ 
ble  eigenvalue  a  for  each  frequency  w.  Many  important  compressible  flows  (e.g., 
high-speed  jets)  have  several  unstable  modes  that  may  be  important  [Sandham 
k  Reynolds,  1989],  and  it  is  desirable  to  be  able  to  characterize  when  multiple 
unstable  inodes  might  arise. 

Unfortunately,  for  the  compressible  Rayleigh  equation,  no  general  results  are 
known — the  number  of  unstable  modes  is  usually  determined  empirically,  by  solv¬ 
ing  the  equations  numerically  and  searching  for  eigenvalues.  However,  for  the 
temporal  stability  analysis  of  the  incompressible  Rayleigh  equation,  some  strong 
results  are  available,  and  it  is  plausible  that  for  the  low  Mach  numbers  we  are 
interested  in,  similar  results  are  true  for  the  compressible  equation.  We  briefly 
summarize  the  most  important  of  these  results  for  the  incompressible  Rayleigh 
equation  (2.12).  Except  where  noted,  these  results  concern  the  tempora/ stability 
analysis,  where  we  fix  a  €  R,  and  search  for  eigenvalues  a;  €  C.  In  particular,  we 
focus  on  unstable  eigenvalues  (for  which. Im(u;)  <  0). 

In  his  landmark  paper,  Rayleigh  [1880]  demonstrated  that  a  necessary  condition 
for  unstable  eigenvalues  is  that  U"  must  change  sign  over  the  range  of  integration — 
i.e..  the  velocity  profile  must  have  an  inflection  point.  Minor  refinements  to  this 
necessary  condition  had  been  made  (see  Drazin  k  Reid  [1981]  for  the  details), 
but  the  next  major  advance  was  not  until  80  years  later,  with  Howard's  circle 
theorem  [Howard,  1961].  Letting  c  =  -uj/a,  this  theorem  states  that  unstable 
eigenvalues  c  must  lie  within  the  semicircle  with  diameter  equal  to  the  range  of 
the  \elocit\  profile  0.  A  few  years  later,  Howard  [1964]  proved  the  following 
powerful  result: 

Theorem  2.1  For  a  class  of  velocity  profiles,  the  number  of  unstable  modes  cannot 
exceed  the  number^  of  inflection  points.  The  class  of  velocity  profiles  U{y)  satisfies 

-  Cs)  for  some  piecewise  continuous  non-negative  function 
K{y)  and  for  some  number  c^. 

The  tanh  velocity  profile  we  consider  falls  into  this  category  (with  =  1/2),  and 
possesses  only  one  inflection  point,  so  for  the  incompressible,  temporal  analysis, 
we  are  guaranteed  only  one  unstable  eigenvalue. 


2.5  Singularity  at  the  critical  layer 
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These  results  were  further  generalized  by  Sattinger  [1967],  who  proved  that 
under  fairly  general  conditions  on  the  velocity  profile,  there  axe  at  most  a  finite 
number  of  complex  eigenvalues  (thus  the  continuous  spectrum  is  real),  and  w^ho 
also  generalized  Howard’s  circle  theorem  to  show  that  all  points  in  the  spectrum 
(not  just  unstable  discrete  eigenvalues)  lie  within  the  circle  with  diameter  equal  to 
the  range  of  C7. 

More  recent  results  for  the  incompressible  Orr-Sommerfeld  equation  (w^hich 
includes  viscous  effects)  were  given  by  Salwen  and  Grosch,  who  showed  that  both 
temporal  and  spatial  problems  have  a  continuous  spectrum  [Grosch  &  Salwen, 
1978],  and  who  gave  a  biorthogonality  relation  for  the  adjoint  problem  [Salwen 
k  Grosch,  1981].  We  discuss  the  biorthogonality  relation  for  the  compressible 
equations  in  section  2.6. 

2.5  Singulcirity  at  the  critical  layer 

It  was  mentioned  in  section  2.1  that  the  governing  equation  (2.11)  has  a  singular 
point  where  U(y)  =  —uja.  Although  this  singularity  lies  on  the  integration  path 
(in  the  complex  y-plane)  only  when  uj  and  a  are  both  real,  we  must  still  be  very 
careful  about  the  treatment  of  this  singularity  when  we  compute  damped  modes 
(with  Im(Q)  >  0),  as  we  describe  below.  The  idea  is  to  construct  a  damped-mode 
solution  that  is  the  analytic  continuation  of  an  unstable-mode  solution,  following 
the  approach  of  Tam  k  Morris  [1980]. 

Suppose  that  uj  +  a[7(yc)  =  9.  but  U\yc)  ^  0,  so  that  yc  is  a  regular  singular 
point  of  equation  (2.11).  (This  locus  of  points  in  the  fluid  where  y  =  yc'is  usually 
called  the  critical  layer  [Drazin  k  Reid,  1981].)  We  first  examine  the  type  of  singu¬ 
larity  of  the  solution  p(y)  at  yc.  A  standard  Frobenius  series  analysis  [Coddington 
k  Levinson,  1955]  shows  that  the  indicial  equation  has  two  solutions  which  differ 
by  an  integer,  and  so  the  solution  p  has  a  logarithmic  singularity  at  yc-  For  the 
velocity  profile  (2.13),  with  cj  >  0  and  Im(a)  <  0  (an  unstable  mode),  the  singu¬ 
larity  yc  lies  in  the  lower  half  y-plane.  Thus,  since  the  solution  is  valid  for  real  y, 
the  branch  cut  of  the  logarithmic  function  must  extend  from  yc  to  infinity  in  the 
lower  half  y-plane,  as  shown  in  figure  2.1. 

If  we  increase  the  frequency  a;,  the  locus  of  yc  moves  toward  the  real  axis  of  the 
y>plane.  and  eventually  crosses  it  (when  Im(Q)  >  0,  thus  becoming  a  stable  mode). 
Once  yc  crosses  the  real  axis  into  the  upper-half  plane,  we  must  be  careful  not  to 
integrate  through  the  branch  cut  which  extends  from  yc  into  the  lower-half  plane. 
For  this  damped-mode  solution  to  be  the  analytic  continuation  of  the  unstable¬ 
mode  solution,  we  must  deform  the  integration  contour  into  the  upper-half  plane 
to  avoid  crossing  the  branch  cut  (see  figure  2.1). 

When  deforming  the  contour  in  the  complex  y-plane,  we  must  pay  close  atten¬ 
tion  to  the  domain  of  the  velocity  profile  C7(y),  to  avoid  inadvertently  integrating 
around  a  branch  point  of  the  map  U  ^  y.  For  the  tanh  profile  (2.13),  the  function 
U  ^  y  Is  multiple-valued,  and  has  two  branch  points,  the  points  0  and  1  on  the 
real  sods,  as  shown  in  figure  2.2.  The  usual  integration  contour  (on  the  real  y-axis) 
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(a)  Unstable  eigenvalue  (b)  Stable  eigenvalue 


Figure  2.1:  Location  of  branch  points  of  the  solution  p{y),  and  integration  contours 
for  the  Rayleigh  equation. 


Figure  2.2:  Domain  of  velocity  profile  U{y),  with  branch  cuts  indicated  in  the 
f/-plane.  Howard’s  circle  theorem  states  that  eigenvalues  —ufa  must  lie  within 
the  dark  shaded  circle  in  the  t7-plane,  which  maps  to  the  dark  shaded  strip  in  the 
t/-plane. 


2.6  Adjoint  equations 
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maps  to  the  interval  (0, 1)  in  the  [/-plane,  so  we  choose  branch  cuts  from  (—oc,  0] 
and  [l,oo)  in  the  [/-plane,  as  shown.  The  range  of  this  branch  of  the  function 
U  y  is  the  strip  -7r/4  <  lm{y/6u})  <  7r/4.  One  must  be  careful  not  to  deform 
the  integration  contour  outside  this  strip,  or  the  corresponding  contour  in  the  U- 
plane  will  go  through  the  branch  cuts,  and  go  the  wrong  way  around  the  branch 
points,  changing  the  value  of  the  integral.  Note  also  that  by  Howard’s  circle  theo¬ 
rem,  eigenvalues  —co/a  (and  therefore  critical  points  where  U{yc)  =  —co/a)  must 
lie  within  the  circle  indicated,  which  maps  to  the  strip  — 7r/8  <  lm{y/S^^)  <  tt/S  in 
the  2/-plane. 

Profiles  from  data.  Sometimes  it  is  desirable  to  use  a  velocity  profile  deter¬ 
mined  from  data,  either  experimental  or  numerical.  In  this  case,  in  order  to  indent 
the  contour  into  the  complex  plane,  one  must  fit  an  analytic  function  to  the  data, 
and  verify  that  the  integration  contour  does  not  leave  the  region  of  convergence  of 
the  function.  In  this  thesis,  whenever  we  have  stable  (or  nearly  stable)  eigenvalues, 
we  fit  our  measured  profiles  to  tanh  profiles,  to  avoid  these  difficulties. 

2.6  Adjoint  equations 

To  determine  how  disturbances  are  amplified  by  the  shear  layer,  we  shall  want 
to  expand  an  arbitrary  disturbance  as  a  sum  of  eigenfunctions  of  the  Rayleigh 
problem.  (It  is  shown  in  Salwen  Grosch  [1981]  that  the  temporal  eigenfunctions 
form  a  complete  set,  but  their  result  does  not  hold  for  the  spatial  eigenfunctions.) 
In  order  to  obtain  the  coefficients  in  this  expansion,  it  is  useful  to  have  a  set  of 
functions  which  are  orthogonal  to  the  eigenfunctions.  We  may  obtain  such  a  set 
of  functions  by  solving  the  adjoint  problem. 

We  begin  defining  the  inner  product  on  complex  vector- valued  functions 

(•.•):L2(r)x/.^(r)-Cby 

{f,9)  =  l^g'fdy.  (2.26) 

where  *  denotes  Hermitian  transpose,  and  T  is  the  domain  of  the  ODE  we  are 
considering.  Here,  we  take  F  =  [-1,  oc),  but  the  extension  to  contours  off  the  real 
axis  is  trivial.  Under  this  inner  product,  the  formal  adjoint  of  is 

L.a 

Writing  =  {u,v,p)^,  the  Lagrange  identity  then  gives 


f —p{iuj* -h  ia*U)  0  — za*  ^ 

pU’  +  m*£7)  -I; 

-m*  -A/2(ic..*  +  ta*t7)y 


(2.27) 
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for  any  tp,  <p  (see  Coddington  &  Levinson  [1955]  or  Stakgold  [2000]  for  a  general 
description  of  the  Lagrange  identity).  Integrating  over  the  domain  F,  we  obtain 


9)  -  =  {v*p  +  ^v) 


-1 


(2.29) 


Note  that  the  boundary  terms  vanish  if  we  impose  the  boundary  conditions 


v(-l)  =  0,  (2.30) 

p^O  asy->oo.  (2.31) 

Gaussian  elimination  on  the  equation  =  0,  again  yields  a  single  equation  for 

the  adjoint  pressure 


dy^‘ 


P  dy 


-j - p'  -  pM^{uj  +  aUf  + 


p'  aU' 
p  u!  +  aU 


+  2 


aU' 


u)  +,  all 


as  well  as  the  expressions 


qU" 
u)  +  aU 
(2.32) 


u  = 


p{u!  +  aU) 

-1 


V  — 


-ry~.  dp* 

.  f-,.  ,  pU  u - — 

p[iijj  +  lab )  \  dy 


Using  (2.34).  the  boundary  condition  (2.30)  becomes 


(2.33) 

(2.34) 


d^  ccU'  ^ 

-j—  - - f-p  at  y  =  -1 

dy  u}  +  aU 


(2.35) 


The  farficid  boundary  condition  is  the  same  as  (2.19),  with  p  replaced  by  p'.  The 
resulting  two-point  BVP  is  solved  as  described  in  section  2.3  for  the  Rayleigh 
equation. 


2.7  Biorthogonality  relation 

We  now  establish  the  following  relation  between  solutions  of  the  Rayleigh  equation 
and  adjoint  solutions. 

Theorem  2.2  (Biorthogonality)  =  0  and  =  0,  then 

(u;i -tj2)  ((91.92))  +  (qi  -  02)  |(pi,ip2l  =  0,  (2.36) 


2. 7  Biortbogonality  relation 
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where  {{f,g))  =  {Af,g)  and  |/,pl  =  {Bf,g),  with 

fp  0  Q\  fpU  0  1  \ 

Aiy)  =  0  p  0  ,  Biy)  =  \  0  pU  0  .  (2.37) 

Vo  0  Af7  V  1  0 

Proof:  We  have 

0—  ,ai  1  *^cj2,Q!2  9^2  ^ 

—  't'u'2,a2 )  9^1  j  *^2) 

=  2(0^1  -UJ2)  {Aipi.^2)  +  i{o:i  -a2)  {B^i,^2) 

=  i{ujl  -UJ2)  ((9^1,  ?2))  -  0^2)  bl,<?2l  , 

which  establishes  the  result.  | 

Strictly  speaking,  the  term  “orthogonal”  is  not  entirely  appropriate  here:  although 
((',•))  is  indeed  an  inner  product  (assuming  p,M  >  0),  on  the  contrary  I*,*!  is 
t>"pically  not  positive  definite  (e.g.,  for  \pU\  <  1  and  |M^C7|  <  1),  and  thus  is  not 
truly  an  inner  product.  However,  we  retain  the  terminology  to  be  consistent  with 
the  literature  [Salwen  &  Grosch,  1981;  Colonius  et  ah,  1998]. 

The  two  biorthogonality  relations  follow  immediately: 

Corollary  2.3  If  =  0  and  =  0;  with  uJi  ^  w)2,  then  ((9:^1,  ^^2))  =  0- 

Corollary  2.4  =  0  and  a\  ^02,  then  |[(^i,92l  —  0- 

Example:  expansion  in  spatial  eigenf^ unctions.  Suppose  we  are  given  a  sinu¬ 
soidal  disturbance,  imposed  at  a  fixed  streain'wise  location  xq.  with  some  known  fre¬ 
quency  and  we  wdsh  to  determine  how  this  disturbance  evolves  as  it  propagates 
downstream.  Specifically,  write  the  disturbance  as  {u,i\p){xo,y,t)  =  /(y)e^“’'^ 
where  f{y)  =  (ti.  t’,p)(xo,  y,  0).  To  solve  this  problem,  we  expand  f{y)  in  terms 
of  (spatial)  eigenfunctions  of  the  Rayleigh  problem.  We  can  then  easily  determine 
how  the  disturbance  evolves,  since  we  know  how  the  eigenfunctions  evolve. 

Assuming  the  spatial  eigenfunctions  form  a  complete  set.  we  expand  f{y)  as  a 
sum  of  eigenfunctions  with  different  values  of  a: 

/(y)  =  XI  +  /  a{i')Pu.-Ay)dn,  (2.38) 

j=l  Ja{uj) 

where  N(uj)  is  the  number  of  discrete  eigenvalues  for  the  given  u;,  cr{uj)  is  the 
continuous  spectrum  of  L^,q.  and 


j  =  1,... 
u  E 


(2.39) 

(2.40) 
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Suppose  Qi  is  the  only  unstable  eigenvalue,  and  we  axe  interested  in  the  coefficient 
ui  of  the  corresponding  eigenfunction,  to  see  how  this  disturbance  will  grow  as  it 
propagates  downstream.  Then,  to  obtain  the  coefficient  of  the  eigenfunction  q,  , 
we  find  the  corresponding  adjomt  eigenfunction  which  satisfies  = 

0.  Then,  applying  theorem  2.2,  we  obtain 


ai  ,  (2.41) 

which  is  easily  solved  for  oi  (provided  the  coefficient  on  the  left-hand  side  is  non¬ 
zero,  which  is  the  generic  case). 


Chapter  3 

Physics  of  Cavity  Oscillations 


This  chapter  presents  a  detailed  study  of  the  physical  mechanisms  underlying 
oscillations  in  the  flow  past  a  cavity.  Direct  numerical  simulations  are  used  to 
study  the  flow,  and  reveal  two  distinct  types  of  oscillations:  a  shear-layer  mode, 
which  resembles  most  cavity  flows  which  have  been  studied  in  the  past,  and  a 
wake  mode,  which  is  less  well  understood.  We  present  a  model  of  the  scaling  laws 
governing  the  oscillations,  and  use  the  model  to  predict  transitions  bet'ween  the 
flow  regimes. 

3.1  Numerical  method 

Flows  with  self-sustained  oscillations  are  inherently  difficult  to  model  computation¬ 
ally,  because  of  their  sensitivity  to  disturbances.  In  the  cavity  flow,  for  instance, 
the  shear  layer  spanning  the  cavity  amplifies  flow  disturbances,  so  small  numeri¬ 
cal  errors  at  the  cavity  leading  edge  will  be  amphfied  exponentially  as  they  travel 
downstream.  If  these  numerical  errors  arise  from  artificial  reflections  at  a  com¬ 
putational  boundary,  for  instance,  they  may  be  indistinguishable  from  physical 
disturbances,  and  may  even  cause  the  flow  to  resonate  at  non-physical  frequencies. 
Furthermore,  in  the  cavity  flow,  the  governing  feedback  mechanism  is  acoustic, 
so  a  high-order,  low-dissipative  numerical  method  is  necessary  to  accurately  re¬ 
solve  acoustic  waves,  which  are  usually  many  orders  of  magnitude  smaller  than 
hydrodynamic  flow  disturbances. 
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3,1.1  Equations  of  motion 

W'e  solve  the  fully  compressible  Navier-Stokes  equations  in  two  dimensions,  written 
in  conservative  form  as  follows: 


dp  d 


dpui 


1  d  f  dui  duj 


—  +  —{puiUj  +p6ij)  =  — — (  _  +  ^  _ 


dx 


dxi 


Re  dxj  \  dxj 

a  -  ^(('  +  1’)“.) '  +  S;  - 


2  duk , 

sax,.  '■ 


together  with  the  equation  of  state 


+ 


1 


Re  Pr  dx^dxi; 

(3.1) 


7-1  rr. 

P  =  ^ - pT. 

7 

Here,  p  is  the  density,  n,-  the  velocity  in  the  direction  of  coordinate  x;,  p  the  pres¬ 
sure.  T  the  temperature,  7  the  ratio  of  specific  heats,  and  e  the  energy,  defined  by 
^  Pi^  —  l^!*^/2)-  w’here  E  is  the  internal  energy  per  unit  mass.  These  equations 
ha\e  been  nondimensionalized  according  to  the  usual  compressible  nondimension- 
alization 


II 

Q 

X,  =  — 

t  = 

Poc 

'  D 

D 

XLi 

n- 

T  = 

T% 

^oc 

r  0 

Poc^oc 

_ 

p  —  - 

—  P^^ocD 

II 

CpP 

where  the  superscript  0“^  denotes  the  dimensional  quantity.  Here.  Cp  is  the  specific 
heat  at  constant  pressure,  k  the  thermal  conductivity,  Re  the  Reynolds  number, 
and  Pr  the  Prandtl  number.  The  length  scale  D  is  the  cavity  depth.  In  this  thesis^ 
we  will  usually  refer  to  the  Revmolds  number  based  on  velocity  Ux  and  momentum 
thickness  6q  at  the  cavity  leading  edge,  defined  by 


P 


M 

L/9o 


Re. 


As  temperature  differences  are  small,  we  assume  constant  transport  properties. 
^^o  take  Pr  —  0./  and  7  —  1.4,  the  values  for  air.  These  equations  are  solved 
directly,  meaning  we  use  no  turbulence  model,  and  resolve  all  of  the  scales  of  the 
flow. 


3.1  Numerical  method 
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3.1.2  Numerical  scheme 

The  equations  (3.1)  are  solved  using  a  sixth-order  compact  finite-difference  scheme, 
described  in  Lele  [1992],  with  a  fourth-order  Runge-Kutta  time  advancement.  This 
combination  of  schemes  results  in  very  low  dispersion  and  numerical  dissipation, 
which  allows  for  accurate  wave  propagation.  The  scheme  has  been  used  previously 
to  study  sound  generation  in  mixing  layers  and  jets,  and  is  capable  of  resolving 
acoustic  fields  with  velocity  fluctuations  five  orders  of  magnitude  smaller  than 
near-field  fluctuations  [Colonius  et  al.,  1997].  The  method  relies  solely  on  physical 
viscosity  for  stabihty.  A  Cartesian  grid  is  used,  with  clustering  of  grid  points  near 
the  walls,  and  in  the  shear  layer  spanning  the  cavity.  Analytical  error  function 
mappings  are  used  for  the  grid  stretching. 

Boundary  conditions  play  a  key  role  in  aeroacoustic  computations.  Artificial 
boundaries  (inflow/outflow/normal)  must  allow  vortical  and  acoustic  waves  to  pass 
freely  with  minimal  reflection.  The  exact  nonreflecting  boundary  conditions  for 
multidimensional  problems  are  nonlocal  in  both  space  and  time  (i.e.,  one  must 
solve  convolution  integrals),  so  in  practice  one  typically  uses  local  approximations, 
which  take  the  form  of  partial  differential  equations  to  be  solved  at  the  boundaries, 
and  errors  in  these  approximations  lead  to  reflections.  Furthermore,  another  type 
of  reflections  arises  from  dispersive  effects  of  the  finite-difference  scheme  used. 
These  errors  are  usually  referred  to  as  “spurious”  numerical  reflections,  and  are 
poorly  resolved  on  the  grid,  while  the  former  reflections  are  well  resolved,  and 
arc  referred  to  as  “smooth”  reflections.  For  linear  equations,  it  is  possible  to 
derive  local  boundary  conditions  that  are  nonreflecting  to  arbitrary  high  order  of 
accuracy,  for  both  types  of  reflections  [Rowley  &:  Colonius,  2000],  but  for  nonlinear 
equations  the  methods  fail.  At  inflow  and  normal  boundaries,  one  might  be  able  to 
linearize  the  boundary  conditions  about  the  mean  flow,  since  the  fluctuations  are 
small,  but  at  the  outflow,  flow  disturbances  are  large,  severely  limiting  the  accuracy 
of  the  linear  boundary  conditions  [Colonius  et  al.,  1993].  For  nonlinear  problems, 
recent  techniques  involve  using  a  robust  but  low-order  boundary  condition  (e.g., 
Thompson  [1987],  Poinsot  &  Lele  [1992]),  in  conjunction  with  a  buffer  zone  near  the 
computational  boundary,  to  attenuate  the  reflections.  These  buffer  zone  methods 
involve  combinations  of  grid  stretching  and  filtering  [Colonius  et  al,  1993],  and 
the  addition  of  forcing  terms  to  the  equations  [Freund,  1997]. 

The  present  method  uses  the  boundary  conditions  of  Poinsot  Lele  [1992] 
for  inflow,  outflow,  and  normal  boundaries,  along  with  the  buffer  zone  method  of 
Freund  [1997].  In  the  buffer  zone,  forcing  terms  of  the  form  a(q  —  qtarget) 
added  to  the  right-hand  side  of  the  equations  (3.1),  where  q  =  {p,pu,pv,e)  is  the 
vector  of  conservative  variables,  and  qtarget  desired  mean  flow.  The  gain  a 

varies  smoothly  from  a  constant  value  at  the  boundary,  to  zero  at  the  edge  of 
the  buffer  (i.e.,  the  edge  of  the  “physical”  portion  of  the  computational  domain). 
The  wall  boundary  conditions  are  isothermal,  at  the  same  temperature  as  the 
freestream.  using  the  formulation  recommended  by  Poinsot  &  Lele  [1992].  Since 
temperature  variations  are  small,  transport  properties  are  assumed  constant,  and 
the  Prandtl  number  is  taken  as  0.7,  its  value  for  air. 
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Figure  3.1:  Schematic  diagram  of  cavity  configuration  and  computational  domain. 

A  schematic  diagram  of  the  computational  diagram  is  shown  in  figure  3.1. 
Tj-pical  grids  (see  captions  of  figures  3.2  and  3.3)  contain  about  half  a  million  grid 
points.  The  code  is  parallelized  using  a  domain  decomposition  method,  and  has 
been  run  on  8  to  32  nodes  of  an  IBM  SP2. 

The  initial  condition  for  all  simulations  is  a  (laminar)  flat-plate  boundary  layer 
along  the  wall,  and  spanning  the  cavity.  We  use  the  (incompressible)  Blasius 
profile  and  spreading  rate  [W’hite,  1991].  The  following  parameters  may  be  varied 
independently:  the  cavity  length,  relative  to  the  initial  momentum  thickness  at  the 
cavity  leading  edge.  L/ffo:  the  Refolds  number  Reg  =  PocUeo/fx:  the  freestream 
Mach  number  M  =  U/a^o'.  and  the  cavity  aspect  ratio  L/D.  Because  of  the 
expense,  only  a  small  portion  of  parameter  space  may  be  investigated,  and  table  3.1 
gives  the  parameters  for  the  runs  performed. 

Note  that  the  quoted  values  of  t^o  correspond  to  the  momentum  thickness  of 
the  boundary  layer  at  the  upstream  edge,  for  the  initial  condition.  For  the  runs  in 
shear-la\er  mode,  the  variation  from  this  initial  \'alue  is  small,  but  for  the  runs  in 
wake  mode,  the  variation  is  substantial.  We  use  the  initial  momentum  thickness 
because  it  is  a  parameter  which  is  well  defined  a  priori. 


3.1.3  Validation 

For  flows  with  self-sustained  oscillations,  it  is  particularly  important  to  verify 
that  results  are  independent  of  the  location  of  the  boundaries,  and  the  boundary 
treatment.  As  mentioned  in  the  previous  section,  if  boundaries  are  not  treated 
properly,  repeated  non-physical  reflections  of  waves  can  lead  to  self-forcing  of  the 
flow ,  in  a  process  that  is  ultimately  indistinguishable  from  physical  instability. 
We  have  run  several  cases  with  variable  boundary  locations  and  grid  spacings,  to 
determine  appropriate  boundary  locations  and  demonstrate  grid  convergence.  The 
results  of  tw’o  of  these  tests  are  shown  in  figure  3.2  and  3.3. 

Figure  3.2  show's  the  normal  velocity  at  a  single  point  in  the  shear  layer  (i  = 


3.1  Numerical  method 
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Run 

L/D 

L/Oo 

M 

Ree 

Mode 

St p  {wake  mode) 

LI 

1 

20.3 

0.6 

73.9 

NO 

- 

0.002 

L2 

2 

52.8 

0.6 

56.8 

SL 

- 

0.008 

L3 

3 

75.0 

0.6 

60.0 

M 

- 

0.031 

L4 

4 

102.1 

0.6 

58.8 

•w 

0.064 

0.227 

L5 

5 

123.3 

0.6 

60.8 

w 

0.054 

0.404 

4M2 

4 

102 

0.2 

58.8 

SL 

- 

4M3 

4 

102 

0.3 

58.8 

M 

- 

4M4 

4 

102 

0.4 

58.8 

W 

0.064 

4M5 

4 

102 

0.5 

58.8 

w 

0.064 

4M6  (L4) 

4 

102 

0.6 

58.8 

w 

0.064 

4M7 

4 

102 

0.7 

58.8 

w 

0.061 

4M8 

4 

102 

0.8 

58.8 

w 

0.061 

2M2 

2 

52.8 

0.2 

56.8 

NO 

- 

2M3 

2 

52.8 

0.3 

56.8 

NO 

- 

2M4 

2 

52.8 

0.4 

56.8 

SL 

- 

2M5 

2 

52.8 

0.5 

56.8 

SL 

- 

2M6  (L2) 

2 

52.8 

0.6 

56.8 

SL 

- 

2M7 

2 

52.8 

0.7 

56.8 

SL 

- 

2M8 

2 

52.8 

0.8 

56.8 

SL 

- 

HI 

1 

23.2 

0.6 

86.3 

NO 

- 

H2 

2 

58.4 

0.6 

68.5 

SL 

- 

H3 

3 

84.9 

0.6 

70.6 

M 

- 

H4 

4 

116.7 

0.6 

6.8.5 

W 

0.063 

LG6a 

6 

45.18 

0.6 

58.57 

NO 

- 

LG6b 

6 

90.36 

0.6 

29.28 

SL 

- 

LG8 

8 

60.24 

0.6 

58.57 

NO 

- 

TK4a 

4 

30.12 

0.6 

58.6 

NO 

- 

TK4b 

4 

60.24 

0.6 

58.8 

SL 

- 

TK4c 

4 

75.30 

0.6 

58.8 

M 

- 

TK4d  (L4) 

4 

102.1 

0.6 

58.8 

W 

0.064 

R4a 

4 

86.06 

0.6 

45.8 

SL 

- 

R4b  (TK4b) 

4 

60.24 

0.6 

58.8 

SL 

- 

R4c 

4 

74.55 

0.6 

80.5 

W 

0.063 

Table  3.1:  Parameters  for  the  different  computer  runs.  Abbreviations  for  modes 
are:  NO  =  No  oscillations,  SL  =  Shear  Layer,  W  =  Wake,  M  ==  Mixed.  Strouhal 
numbers  of  vortex  shedding  (based  on  cavity  depth  and  freestream  velocity)  are 
given  for  wake  mode  runs,  and  the  computed  drag  coefficient  is  given  for  run  series 
L  (see  section  3.3). 
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Figure  3.2.  Effect  of  boundary  position  and  grid  resolution  on  the  normal  velocity 

at  y  0  and  x  =  3.13Z).  (a)  shows  reference  case  L4  ( - )  compared  with 

a  finer  grid  case  ( - )  and  larger  domain  case  (• . );  (b)  shows  only  the  larger 

domain  case,  for  longer  time.  Reference  case  L4  has  downstream  boundary:  10. 6D, 
upstream:  -4.3£>,  normal:  9.2D.  Grid  has  1152  x  384  points  above  the  cavity  in 
the  streamwise  and  spanwise  directions,  respectively,  and  384  x  94  points  in  the 
cavity.  Finer  grid  case  has  same  boundaries  as  run  L4,  but  50%  more  grid  points 
(in  each  direction).  Larger  domain  case  extends  to  15D  downstream  and  15D  in 
the  normal  direction.  Note  that  the  dotted  line  falls  nearly  directly  on  top  of  the 
solid  line  in  (a). 


3.13£).  y  —  0).  for  different  grid  resolutions  and  boundary  locations,  but  otherwise 
using  the  same  parameters  as  run  L4.  Other  probe  locations  yielded  similar  results. 
Note  that  time  is  normalized  by  the  freestream  velocity  and  the  total  length  of 
the  computational  domain  for  the  reference  case.  The  larger-domain  case  is  run 
for  longer  than  the  finer-grid  case,  as  the  domain  size  is  more  hkely  to  cause  low 
frequency  errors,  arising  from  nonphysical  reflections  at  boundaries.  Figure  3.3 
shows  a  similar  test  for  run  L2,  and  the  spectra  are  compared. 

Over  3  to  4  flow-through  times,  the  results  are  nearly  identical,  independent 
of  grid  resolution  and  boundary  location.  Small  differences  are  apparent  at  later 
times,  which  is  not  unexpected  given  the  sensitive  dependence  on  initial  conditions. 
For  run  L2.  greater  differences  are  observed  after  about  10  flow-through  times,  but 
the  spectra  are  almost  identical.  We  conclude  that  the  boundary  locations  and  grid 
resolutions  for  runs  L2  and  L4  are  adequate,  and  similar  locations  and  resolutions 
were  used  in  the  other  runs  in  table  3.1. 

3.2  Shear-layer  mode 

The  shear-layer  mode  is  characterized  by  the  feedback  process  described  in  the 
introduction:  the  roll-up  of  vorticity  in  the  shear  layer,  impingement  and  scat¬ 
tering  of  acoustic  waves  at  the  downstream  cavity  edge,  upstream  acoustic  wave 
propagation,  and  receptivity  of  the  shear  layer  to  acoustic  disturbances.  The 
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(a)  (b) 


Figure  3.3:  Effect  of  boundary  placement  on  the  normal  velocity  at  y  —  0  and 

X  =  1.57Z),  for  run  L2.  (a)  shows  the  time  trace  for  the  reference  case  ( - ) 

and  for  a  larger  domain  case  ( . );  (b)  shows  the  spectra  of  the  data  in  (a). 

Reference  case  L2  has  downstream  boundary  7.6D,  upstream:  ~-3.9D,  normal: 
9.2jD.  Grid  has  1008  x  384  gridpoints  above  the  cavity,  and  240  x  96  gridpoints 
in  the  cavity.  Larger  domain  case  extends  to  11. 8D  downstream  and  15. 6D  in  the 
normal  direction. 

process  is  clearly  born  out  by  the  computational  results.  Figure  3.4  shows  vor- 
ticity  contours  at  three  different  times  for  run  L2,  and  these  are  indicative  of  all 
the  runs  in  shear-layer  mode  of  oscillation.  Vortical  disturbances  propagate  and 
grow  along  the  shear  layer,  while  the  flow  inside  the  cavity  is  relatively  quiescent. 
A  weak,  relatively  steady  vortex  occupies  the  downstream  half  of  the  cavity,  as 
found  by  Roshko  [1955].  The  steady  vortex  induces  vorticity  of  the  opposite  sign 
(to  boundary-laver  vorticity)  dong  the  cavity  walls.  The  steadiness  of  the  vortex 
suggests  that  the  interaction  of  the  flow  inside  the  cavity  with  the  shear  layer  is 
relatively  weak.  The  shear  layer  will  be  examined  in  more  detail  in  sections  3.2.3 
and  3.2.4. 

3.2.1  Acoustic  field 

In  general,  it  is  not  possible  to  distinguish  objectively  between  acoustic  and  hydro- 
dynamic  fluctuations  in  a  complex,  vortical  flow.  In  this  section,  we  examine  the 
acoustic  field  produced  by  the  cavity,  where  we  loosely  define  the  acoustic  field  as 
the  irrotational  field  generated  as  vorticity  is  swept  past  the  trailing  corner  of  the 
cavity.  For  ver\'  low  Mach  numbers,  these  disturbances  are  essentially  incompress¬ 
ible,  and  the  propagation  speed  becomes  arbitrarily  large  as  the  Mach  number 
goes  to  zero.  Here  we  consider  Mach  numbers  from  0.2  to  0.8,  and  simply  refer  to 
the  scattered  field  as  ‘‘acoustic,'’ 

The  density  fluctuations  produced  by  the  self-sustained  oscillations  is  shown 
in  figure  3.5,  where  we  compare  the  present  results  to  schlieren  photographs  taken 
by  Krishnamurty  [1956]  for  Mach  numbers  0.64,  0.7,  and  0.8,  and  LjD  —  2.  The 
experimental  conditions  are  similar,  in  that  the  incoming  boundary  layer  is  lami- 
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Figure  3.4.  Instantaneous  vorticity  contours  for  run  L2  (shear-layer  mode)  at  3 
different  times  (a-c),  corresponding  to  approximately  one-third  phase  intervals  of 
the  dominant  mode  of  oscillation  (2nd  Rossiter  mode).  15  equi-spaced  contours 
between  ujD/U  =  -5  and  1.67  are  shown;  positive  contours  are  dashed.  Only  a 
small  portion  of  the  computational  domain  near  the  cavity  is  shown. 


nar,  but  the  Reynolds  number  in  the  experiment  is  higher  by  a  factor  of  about  5. 
The  cavity  in  the  experiment  had  a  transverse  width  of  4  inches,  compared  to  the 
depth  of  0.1  inch,  and  it  is  expected  under  these  conditions  that  the  instabilities 
are  approximately  two-dimensional.  There  is  a  very  good  qualitative  agreement 
between  experiment  and  computation.  The  greyscale  contour  levels  in  the  simu¬ 
lation  where  chosen  to  provide  the  best  visual  comparison  with  the  experiment. 
Note  from  the  figure  that  the  frequencies  of  oscillation  (and  hence  wavelength  of 
the  acoustic  field)  are  quite  close  at  M  =  0.7  and  M  =  0.8.  At  M  0.6  (0.64  for 
the  experiment),  it  appears  that  the  frequency  of  oscillation  is  substantially  low'er 
in  the  simulation,  and  the  experiment  at  these  conditions  is  dominated  by  Rossiter 
mode  2  while  the  computation  is  dominated  by  Rossiter  mode  1  (see  section  3.2.2). 

Also  note  from  figure  3.5  that  the  radiation  is  intense  enough  that  there  is 
significant  steepening  of  the  compressions  as  the  w'av'es  propagate,  especially  at 
the  higher  Mach  numbers. 

Figure  3.6  shows  the  sound  pressure  levels  (SPL)  for  the  acoustic  field  above 
the  cavity,  for  A/  =  0.6.  The  maximum  SPL  is  very  high,  about  180  dB,  at  a 
point  near  the  dow-nstream  edge,  and  peak  radiation  to  the  far  field  occurs  at 
an  angle  of  about  135  degrees  from  the  flow  direction.  In  this  direction,  sound 
pressure  levels  still  reach  170 dB  at  a  distance  of  3  cavity  depths.  These  numbers 
are  quite  high,  but  Krishnamurty  [1956]  estimated  sound  pressure  levels  for  the 
case  of  laminar  boundary  layers  in  excess  of  163  dB  for  a  variety  of  geometries. 
These  estimates  were  based  on  deflections  from  finite-fringe  interferometry,  but 
no  detailed  mapping  to  SPL  was  performed.  Krishnamurty  [1956]  found  that 
flows  with  laminar  boundary  layers  produce  higher  SPLs  than  flows  with  turbulent 
boundary  layers,  and  that  the  L/D  =  2  cavity  was  louder  than  longer  cavities. 

Figure  3.7  illustrates  the  separate  contributions  to  the  acoustic  field,  when  two 
resonant  frequencies  are  present  simultaneously.  Density  fluctuations  are  plotted 
for  run  L2.  w'hich  has  two  strong  resonant  frequencies,  at  St  =  fL/U  =  0.4  and  0.7. 
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(a)  Schlieren,  A/  =  0.64  (b)  Run  2M6,  M  —  0.6 


(c)  Schlieren,  A/  =  0.7  (d)  Run  2M7,  M  =  0.7 


(e)  Schlieren.  A/  =  0.8  (f)  Run  2M8,  A/  =  0.8 


Figure  3.5:  Comparison  of  schlieren  photographs  [Krishnamurty,  1956]  with  con¬ 
tours  of  density  gradient  from  the  DNS.  In  the  schlieren  photographs,  the  knife 
edge  is  horizontal  in  (c),  vertical  in  the  others;  in  the  DNS  figures,  dpfdy  is  shown 
in  (d),  dpjdx  in  the  others. 
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X 


Figiire  3.6:  Sound  pressure  level  (SPL)  for  run  L2. 

\Vc  performed  a  Discrete  Fourier  Transform  (DFT)  of  125  samples  (every  250  time 
units)  of  the  computational  data  over  a  period  of  time  TU/L  =  30,  corresponding 
to  12  periods  of  the  lower  frequency,  and  21  periods  of  the  higher  frequency.  The 
resulting  data  record  is  approximately  periodic  in  time,  and  the  data  is  detrended 
prior  to  taking  the  DFT.  We  experimented  with  different  windowing  techniques, 
signal  durations,  and  sampling  rates,  and  we  believe  the  results  presented  here  are 
free  from  any  artifacts  of  the  hmited  signal  duration  and  relatively  low  sampling 
rate. 

Real  and  imaginary  parts  of  the  DFT  are  plotted,  and  the  different  wavelengths 
of  farfield  acoustic  radiation  are  apparent.  Hydrodynamic  disturbances  in  the  shear 
layer  are  also  evident,  and  will  be  discussed  further  in  section  3.2.4.  It  is  striking 
that  the  level  of  the  density  fluctuations  in  the  farfield  is  comparable  to  the  level 
of  the  density  fluctuations  in  the  shear  layer.  This  indicates  the  high  efficiency 
of  the  acoiLstic  scattering  process  at  the  cavity  trailing  edge.  The  wavelength  of 
the  acoustic  fields  at  differing  angles  agrees  with  the  predicted  value  based  on  the 
frequency  of  oscillation  and  the  speed  of  the  wavefront,  aoc(l  +  M cos  9),  where  6 
is  the  angle  of  the  wave  propagation  direction  measured  fi'om  downstream  i-axis, 
and  Occ  is  the  ambient  sound  speed. 

3.2.2  Frequencies  of  oscillation 

Figure  3.8  shows  the  spectrum  of  normal  velocity  at  a  point  along  the  shear  layer, 
at  y  =  0,  i/L  =  0.783,  for  a  L/D  =  2  cavity  with  M  =  0.7  (run  2M7),  which 
is  oscillating  in  shear-layer  mode.  Several  distinct  peaks  and  their  harmonics 
are  evident.  The  frequencies  labeled  “Mode  I,  11”  correspond  to  the  first  two 
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Figure  3.8:  Spectrum  of  normal  velocity  at  y  =  0  and  x  =  1.57D,  for  run  2M7, 


frequencies  predicted  by  Rossiter’s  semi-empirical  formula  ' 


_  f„L  _  n  -  7 
U  ~M  +  1/k' 


n  =  l,2,... 


(3.2) 


where  Stn  is  the  Strouhal  number  corresponding  to  the  nth  mode  frequency  /„. 
and  K  and  7  are  empirical  constants  corresponding  to  the  average  convection  speed 
of  disturbances  in  the  shear  layer,  and  a  phase  delajt  The  other  peaks  in  figure  3.8 
are  either  harmonics  or  nonlinear  interactions  between  the  Rossiter  modes. 

In  figure  3.9.  the  frequencies  of  the  two  most  energetic  peaks  in  the  spectra 
for  the  series  of  runs  2iM  with  L/D  =  2  and  L/Oq  =  52.8  and  4M  with  L/D  =  A 
and  1/00  =  102  are  compared  to  experimental  data  and  predictions  from  equa¬ 
tion  (3.2).  ■wnh  7  =  0.25  and  k  =  1/1.75,  the  original  values  used  by  Rossiter. 
The  conditions  in  Krishnamurty’s  1956  experiments  are  closest  to  those  in  our 
simulations,  in  that  the  upstream  boundary  layers  were  laminar,  but  the  Rejmolds 
number  is  fi\e  times  larger  than  that  in  our  simulations.  Krishnamurtv  found  only 
Mode  II  oscillations,  and  the  frequencies  are  slightly  higher  than  those  measured 
in  our  simulations.  Sarohia  [1975]  found  in  his  experiments  in  water  that  the  fre¬ 
quency  increases  with  Reynolds  number.  He  observed  both  Modes  I  and  II,  as 
shown  in  the  figure,  and  it  appears  that  the  frequencies  for  higher  Mach  numbers 
would  agree  well  with  Sarohia’s,  when  extrapolated  to  M  =  0.  For  turbulent 
boundary  layers,  the  frequencies  are  lower  than  for  laminar  boundary  layers,  as 
shown  in  the  figure,  and  as  noted  previously  by  both  Krishnamurty  and  Sarohia. 

Note  that  there  is  quite  a  bit  of  scatter  in  the  experimental  data  at  low  Mach 
numbers.  In  general,  the  data  from  experiments  in  water  are  fairly  consistent 
[Sarohia,  1975;  Knisely  &  Rockwell,  1982],  but  experiments  in  air  at  low  Mach 
numbers  have  more  scatter  [Tam  L  Block,  1978;  Rockwell  k  Schachenmann,  1982; 
Ahuja  i:  Mendoza,  1995].  because  the  frequencies  of  the  resonant  acoustic  modes 
of  the  cavity  become  comparable  to  the  Rossiter  frequencies. 

Peaks  in  the  spectra  corresponding  to  harmonics  and  sum  and  difference  modes 
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X  Karamchetti  (1956) 


Figure  3.9:  Strouhal  numbers  for  peaks  in  spectra  for  the  2M  series  of  runs  (shear- 
layer  mode),  and  4M  series  (shear- layer  and  wake  modes)  compared  to  experiment 
and  equation  (3.2). 


(primarily  harmonics  of  n  =  1  and  the  sum  of  Rossiter  modes  n  =  1  and  n  =  2) 
were  also  detected  in  the  data.  Such  nonlinear  interactions  between  modes  w’ere 
studied  by  Knisely  &:  Rockwell  [1982]  for  incompressible  cavity  flow,  and  more 
recently  by  Cattafesta  et  al.  [1998]  for  moderate  subsonic  speeds,  with  turbulent 
upstream  boundary  layers.  These  studies  used  bicoherence  analysis  and  other 
sophisticated  signal  processing  tools  that  rely  on  long  time  series  of  data.  For 
the  spectra  obtained  here,  signal  durations  were  relatively  short,  on  the  order  of 
20-100  times  the  lowest  frequency  peak  in  the  data,  so  it  was  difficult  to  detect 
the  presence  of  any  mode  switching,  as  was  found  to  occur  over  long  times  by 
Cattafesta  et  al.  [1998].  In  the  computations,  it  appears  that  the  two  primary 
Rossiter  modes  (n  =  1,2)  were  coexistent. 

Finally,  we  note  that*  in  figure  3.9,  oscillation  frequencies  from  the  4M  series  of 
modes  show  the  transition  to  wake  mode  oscillations  for  Al  >  0.3.  These  data  are 
discussed  more  fully  in  section  3.3. 

3.2.3  Shear-layer  spreading  rate 

The  shear  layer  spanning  a  cavity  differs  from  a  mixing  layer  in  two  main  respects: 
Kelvin-Helmholtz  instabilities  are  constantly  being  excited  by  the  intense  acous¬ 
tic  environment,  and  the  entrainment  is  modified  by  the  presence  of  the  cavity. 
Despite  these  differences,  most  researchers  report  that  shear  layers  over  cavities 
closely  resemble  turbulent  free  shear  layers,  in  that  the  spreading  rate  is  approx¬ 
imately  linear.  However,  there  is  some  discrepancy  as  to  the  actual  value  of  the 
spreading  rate. 

Here,  we  use  the  vorticity  thickness  defined  in  section  2.2,  as  a  measure  of 


32 


3  Physics  of  Cavity  Oscillations 


Figure  3.10;  Vorticity  thickness  along  the  shear  layer  for  runs  L2  (a)  and 
LG6b  (□),  and  linear  fits  with  slope  db^jdx  -  0.05  ( - )  and  0.07  ( _ ). 

the  shear-layer  thickness.  It  is  common  to  use  the  momentum  thickness  Q  (also 
defined  in  sec  ‘1.1),  but  the  momentum  thickness  shows  significant  variation  caused 
by  the  recirculating  region  in  the  downstream  portion  of  the  cavity.  The  vortic¬ 
ity  thickness  is  a  local  measure  of  the  maximum  shear,  which  better  determines 
imst ability  properties  of  the  shear  layer,  and  so  it  is  the  thickness  we  choose. 

Sarohia  [19/5]  appears  to  have  been  the  first  to  measure  the  spreading  rate  in 
detail,  and  found  that  the  spreading  was  approximately  linear,  and  the  spreading 
rate  increased  as  L/9q  increased.  As  L/Bo  increased  from  52.5  to  105.2,  with  Rce 
and  D/Oq  held  constant,  the  spreading  rate  db^jdx  varied  from  0.025  to  0.088. 
These  \alues  are  significantly  lower  than  typical  values  for  turbulent  free  shear 
layers,  generally  accepted  to  be  around  db^jdx  0.162  [Brown  fv.  Roshko,  1974). 
The  upstream  boundarj'  layers  in  Harohia’s  experiments  were  laminar.  Gharib 
Aj  Roshko  [1987]  also  found  linear  growth  of  the  shear  layer,  and  found  that  the 
spreading  rate  was  fairly  constant  at  db^^jdx  —  0.124'  for  cavities  with  >  103. 
This  indicates  a  spreading  rate  much  closer  to  that  of  a  turbulent  free  shear  layer. 
The  recent  low  Mach  number  experiments  of  Cattafesta  et  al.  [1997]  also  exhibit 
linear  .spreading,  with  turbulent  upstream  boundary  layers.  Though  the  spreading 
rate  is  not  stated,  they  report  that  for  a  cavity  with  I/^o  =  328,  the  spreading 
rate  closely  matches  that  of  a  turbulent  free  shear  layer,  but  for  a  cavity  with 
Z,/^o  ===  81,  the  spreading  rate  is  507c  higher,  the  opposite  trend  of  that  observed 
by  Sarohia. 

Shear-layer  thicknesses  from  two  of  our  runs  are  plotted  in  figiue  3.10.  Our 
data  also  indicate  approximately  linear  growth,  with  spreading  rates  similar  to 
those  measured  by  Sarohia,  Our  run  L2  (L/^0  =  53)  has  a  spreading  rate  of  about 
db^ldx  =  0.05,  while  Sarohia’s  experiments  had  db^jdx  =  0.056  for  L/^o  =  60. 
Our  spreading  rates  also  increase  with  L/^o-  and  for  run  LG6b  (L/0o  =  90)  we 

‘Where  thicknesses  are  given  in  terms  of  momentum  thickness  6,  we  set  <5,^  =  40,  which  is 
exact  for  a  hyperbolic  tangent  profile. 
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Figure  3.11:  Mean  streamwise  velocity  profile,  u{x^y)/U,  for  run  L2  ( - )  com¬ 
pared  to  the  h\’perbolic  tangent  profiles  with  the  same  vorticity  thickness  ( - ). 

measure  dS^jdx  —  0.07,  compared  to  an  experiment  of  Sarohia’s  with  =  85 
and  dS^fdx  =  0.064.  Thus,  our  measurements  support  the  trend  noted  by  Sarohia, 
that  the  shear  layer  spreads  faster  for  longer  cavities.  A  likely  reason  for  this  trend 
is  that  a  longer  shear  layer  amplifies  disturbances  more,  so  the  final  amplitude  of 
oscillations  is  larger,  thus  increasing  Reynolds  stresses,  which  cause  the  spreading. 

3.2.4  Convection  and  amplification  by  the  shear  layer 

In  this  section,  we  examine  how  disturbances  are  convected  and  amplified  by  the 
shear  layer,  and  we  compare  the  DNS  results  Jo  those  predicted  by  linear  stability 
theory.  The  linear  stability  calculations  are  performed  as  described  in  chapter  2. 
W’e  note  that  previous  analyses  (e.g.,  Tam  k  Block  [1978])  did  not  account  for  the 
cavity  bottom,  as  the  cavity  was  assumed  infinitely  deep  compared  to  the  shear 
layer  thickness. 

W  e  use  the  simulations  to  determine  the  mean  flow  for  the  linear  stability  cal¬ 
culations.  Two  different  calculations  are  performed,  one  using  the  actual  velocity 
profiles  from  the  DNS,  and  another  using  hyperbolic  tangent  profiles  with  the 
same  vorticity  thickness  and  deflection.  As  shown  in  figure  3.11,  the  agreement  is 
good  close  to  the  cavity  leading  edge,  but  much  worse  near  the  rear  of  the  cavity, 
where  the  steady  captive  vortex  is  present.  The  actual  DNS  velocity  profiles  will  of 
course  be  more  accurate,  but  the  tanh  profiles  yield  very  similar  results,  and  later 
(section  3.4)  we  use  the  tanh-profile  linear  stability  calculations  as  a  predictive 
tool  to  predict  the  total  amplification  by  the  shear  layer. 

The  overall  magnitude  and  phase  of  the  instability  wave  is  found  by  integrating 
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the  complex  grow4;h  rate  a  over  the  streamwise  positions; 

f(^<y)  =  fiy)^xp^i  J  a  dx  ),  (3.3) 

where  both  the  wavenumber  a  and  the  eigenfunction  /  eire  slowly  varying  functions 
of  X.  Without  carrying  the  analysis  to  higher  order,  the  normalization  of  the 
eigenfunctions  is  arbitrary,  and  we  have  set  the  maximum  amplitude  of  the  normal 
velocity  mode  to  unit}'  at  each  axial  position. 

Figure  3.12  shows  the  normal  velocity  instability  waves  v,  at  the  two  resonant 
frequencies  present  in  run  L2  (St  =  0.4  and  0.7),  using  the  mean  flow  from  the 
same  run.  Also  shown  is  the  DFT  of  the  normal  velocity,  computed  as  described 
in  section  3.2.1. 

The  linear  stability  eigenfunctions  are  seen  to  be  in  very  good  qualitative  agree¬ 
ment  with  the  Kelvin- Helmholtz  instability  waves  at  the  corresponding  frequencies, 
except  ver}  near  the  cavity  trailing  edge,  where  the  DNS  results  show  some  of  the 
disturbances  being  swept  down  into  the  cavity  by  the  (nearly  steady)  vortex  that 
occupies  the  downstream  half  of  the  cavity. 

A  more  quantitative  measurement  of  the  phase  variation  is  presented  in  fig¬ 
ure  3.13.  The  phase  of  vorticity  disturbances  in  the  shear  layer  is  plotted,  again 
for  run  L2,  at  the  two  resonant  frequencies.  (This  is  just  the  phase  of  the  DFT 
of  vorticity.  along  y  =  0.)  Also  plotted  is  the  phase  of  the  dilatation  close  to 
the  cavity  floor,  which  represents  the  upstream-traveling  acoustic  wave  inside  the 
cavity.  These  variables  were  chosen  to  better  separate  hydrod}'namic  disturbances 
from  acoustic  disturbances,  though  very  similar  results  were  obtained  when  normal 
velocity  and  pressure  were  used  in  the  shear  layer  and  cavity  floor,  respectivelv. 
The  dilatation  phase  is  relatively  constant  in  the  y-direction,  except  near  the  shear 
layer,  where  the  hydrodynamic  fluctuations  are  significant.  Notice  that  the  total 
phase  variation  from  the  shear-layer  convection  (downstream)  and  acoustic  prop¬ 
agation  (upstream)  is  almost  exactly  27rn,  where  n  is  the  index  of  the  Rossiter 
mode.  This  phase  criterion  is  similar  to  that  found  in  several  experiments  [Knisely 
k  Rockwell.  1982:  Rockwell  k  Schachenmann,  1982;  Gharib  k  Roshko,  1987], 
which  show  that  in  the  low  Mach  number  limit,  the  total  phase  variation  in  the 
shear  layer  alone  is  a  multiple  of  2t:.  In  this  limit,  the  acoustic  propagation  is  of 
course  instantaneous,  so  for  our  compressible  simulations  we  must  add  the  phase 
variation  of  the  finite-speed  acoustic  propagation.  Note  that  the  phase  speed  of 
the  acoustic  waves  along  the  cavity  bottom  is  far  from  constant.  This  indicates 
the  presence  of  multiple  acoustic  reflections  by  the  cavity  walls. 

Linear  stability  predictions  for  the  phase  variation  in  the  shear  layer  are  also 
plotted,  and  the  phase  variation  is  well  predicted  except  for  very  near  the  lead¬ 
ing  and  trailing  edges  of  the  cavity,  where  the  DNS  data  show  a  slower  phase 
speed  (steeper  phase  variation).  We  do  not  expect  the  linear  stability  calculations 
to  be  very  accurate  in  these  regions,  because  the  flow  is  rapidly  changing  in  the 
streamwise  direction,  and  significantly  non-parallel  at  the  downstream  corner.  In 
addition,  any  flow/acoustic  coupling  will  be  most  important  near  the  cavity  cor- 
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(d)  DFT  of  DNS,  St2  =  (e)  LS,  DNS  velocity  pro-  (f)  LS,  tanh  velocity  pro- 

0.7  file.  St2  =  0.7  file,  St2  =  0.7 

Figure  3.12:  Comparison  of  mode  shapes  (real  part)  for  normal  velocity  fluctua¬ 
tions  in  the  shear-layer  region  for  parameters  of  run  L2.  Contours  are  between 
±0.01U  for  (a-c),  between  ±0.0051'  for  (d-f). 
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Figure  3.13:  Phase  of  vortical  and  acoustic  disturbances  for  run  L2,  at  the  fre¬ 
quencies  of  first  two  Rossiter  modes  (St  =  0.4,  0.7).  Vorticity  disturbances  in  the 
shear  layer  at  y  =  0  (□,<));  dilatation  along  cavity  bottom  at  y  =  —0.99D  (a,V); 

and  shear-layer  phase  predicted  by  linear  stability  ( - , - ).  Dilatation  phase 

is  shifted  to  line  up  with  vorticity  phase  at  x  =  L. 

ners,  where  the  receptivity  to  acoustic  disturbances  is  high.  These  effects  have 
been  studied  in  detail  by  Rockwell  k  Schachenmann  [1982]. 

Recall  that  Rossiter’s  formula  for  the  resonant  frequencies  includes  an  empirical 
constant  (-y  in  equation  (3.2))  that  represents  an  additional  phase  lag  somewhere 
in  the  feedback  loop.  It  is  possible  that  the  steeper  phase  variation  exhibited  near 
the  leading  and  trailing  edges  of  the  cavity  is  the  cause  of  this  additional  phase 
variation.  The  average  phase  speed  from  the  DNS  data  is  Cp/U  =  0.49  for  St  =  0.4 
and  Cp/L-  =  0.41  for  St  =  0.7,  while  the  linear  stability  calculations  predict  phase 
speeds  of  0.63  and  0.49  for  St  =  0.4, 0.7. 

Figure  3.14  shows  the  amplitude  of  normal  velocity  disturbances  in  the  shear 
layer.  The  measured  growth  rates  are  significantly  smaller  than  those  predicted 
by  the  linear  theory,  which  is  surprising,  because  several  experiments  show  the 
amplitude  is  well  predicted,  at  least  for  moderate  values  of  x/d.  Knisely  k  Rock¬ 
well  [1982]  used  a  constant-thickness  mean  profile,  and  found  that  the  amplitude 
matched  the  linear  theory  well,  for  x/Oq  <  30;  Cattafesta  et  al.  [1997]  found  good 
agreement  for  x/Oq  <  60,  also  using  a  constant-thickness  mean  profile.  However, 
otir  Re\Tiolds  number  is  much  smaller  than  that  in  either  of  these  experiments,  so 
presumably  a  viscous  stability  calculation  would  agree  better. 

In  summary,  linear  stability  theory  gives  reasonable  predictions  for  the  mode 
shapes  of  the  resonant  frequencies,  and  also  the  convection  speeds  of  disturbances, 
but  amplification  rates  are  significantly  over-predicted.  The  linear  stability  calcu¬ 
lation  was  compressible,  but  inviscid,  and  locally  parallel.  Adding  viscous  effects. 
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Figure  3.14:  Amplitude  of  normal  velocity  fluctuations  (run  L2)  along  =  0  at 
Sti  =  0.4  (□)  and  St2  =  0.7  (0);  and  predictions  from  linear  stability  for  Sti 
( - )  and  St2  ( - ). 

including  effects  of  flow/acoustic  coupling,  or  carrying  out  a  multiple  scales  analy¬ 
sis  to  account  for  slightly  non-parallel  effects  [e.g.,  Crighton  k  Gaster,  1976]  may 
provide  a  better  agreement  than  is  obtained  here. 


3.3  Wake  mode 

As  the  length  or  depth  of  the  cavity  (relative  to  the  upstream  boundary-layer 
thickness)  and/or  Mach  and  Reynolds  numbers  are  increased,  there  is  a  substantial 
change  in  the  behavior  of  the  cavity  oscillations.  Under  these  conditions,  the  flow 
is  characterized  by  a  large  scale  shedding  from  the  cavity  leading  edge.  As  noted 
in  the  introduction,  Gharib  k  Roshko  [1987]  were  the  first  to  understand  this 
transition  in  detail,  and  used  the  term  wake  mode  to  describe  the  resulting  flow 
regime.  Connections  with  the  experiment  are  discussed  further  below. 

3.3.1  Flow  features 

Figure  3.15  shows  the  general  features  of  wake-mode  oscillations.  A  vortex  forms 
at  the  leading  of  the  cavity,  and  grows  until  it  is  nearly  the  size  of  the  cavity.  As 
it  is  forming,  irrotational  free  stream  fluid  is  directed  into  the  cavity,  impinging 
on  the  cavity  walls.  Vorticity  of  the  opposite  sign  (to  that  in  the  vortex  and  the 
upstream  boundary  layer)  is  generated  between  the  large  vortex  and  the  upstream 
cavity  wall.  This  opposite  sign  vorticity  accumulates,  and  the  large  vortex  is  shed 
from  the  leading  edge  and  ejected  from  the  cavity  in  a  violent  event.  The  vortex  is 
large  enough  to  cause  flow  separation  upstream  of  the  cavity  during  its  formation, 
and  again  in  the  boundary  layer  downstream  of  the  cavity  as  it  convects  away. 

Time  traces  of  the  normal  velocity  at  a  point  in  the  shear  layer  (y  =  0,  x  = 
3.131?)  are  shown  in  figure  3.16,  for  the  series  of  runs  L1-L5,  where  L/6o  was  varied, 
with  constant  DfOo.  For  the  shortest  cavity  (run  Ll,  L/Oq  =  20),  the  flow  is  steady. 
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Figure  3.15:  Instantaneous  vorticity  contours  for  run  L4  (wake  mode)  at  four 
different  times  (a-d),  corresponding  to  approximately  quarter  phase  intervals  of 
the  periodic  cycle.  15  contours  between  uB/U  =  -5  and  1.67.  Positive  contours 
are  dashed.  Only  a  small  portion  of  the  computational  domain  near  the  cavity  is 
shown. 

For  a  longer  cavity  (run  L2,  L/^o  =  53),  shear-layer  mode  oscillations  are  evident, 
and  as  the  cavity  length  increases  further  (runs  L4,  L5,  L/0o  >  100),  the  flow 
transitions  into  wake  mode  after  a  few  oscillations  in  shear-layer  mode.  For  run  L3 
{L/Oq  75)  it  appears  that  there  is  mode  switching,  with  wake  and  shear-layer 
modes  being  present  at  different  times  (what  we  referred  to  in  table  3.1  as  mixed 
mode).  The  transition  is  also  a  function  of  Mach  number,  and  for  1/00  =  102, 
shear  layer  mode  occurs  for  M  <  0.3,  and  wake  mode  for  M  >  0.3.  Time  traces  in 
the  different  regimes  look  similar  to  those  in  figure  3.16,  and  again  indicated  the 
presence  of  a  mixed  mode. 

Figure  3.17  contrasts  the  time-averaged  flow  for  runs  L2  (shear-layer  mode)  and 
L4  (wake  mode).  The  mean  streamlines  in  wake  mode  are  significantly  deflected 
above  the  cavity,  and  show  that  on  average,  the  boundary  layer  upstream  of  the 
cavity  sees  an  adverse  pressure  gradient.  On  average  the  flow  in  the  cavity  is 
strongl}’  recirculating,  and  there  is  an  impingment  of  the  recirculating  flow  on 
the  rear  wall.  It  is  important  to  contrast  the  mean  flow  with  the  instantaneous 
visualizations  of  figure  3.15,  which  shows  that,  unhke  in  shear-layer  mode,  there  is 
no  stationar3^  vortex  within  the  cavity  instantaneously.  The  region  of  high  pressure 
near  the  back  corner  of  the  cavity  resembles  that  observed  by  Fox  [1965]  in  his 
high-drag  flow  regime.  Variations  in  the  average  coefficient  of  pressure  are  also 
quite  large,  reaching  a  minimum  of  about  -0.5  where  the  flow  is  expanding  into 
the  cavity,  to  about  0.3  in  the  impingment  region  on  the  rear  step. 


3.3  WaJce  mode 


Figure  3.17:  Time-averaged  flow  for  (a)  L2  (shear-layer  mode),  and  (b)  4M6  (wake 
mode).  Mean  streamlines  (solid  lines)  are  superposed  on  contours  of  constant  Cp 
(dashed  lines). 
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Figure  3.18:  Mean  streamwuse  velocity  profiles,  u{x,y)/U,  at  different  streamwise 
positions.  Run  L2  ( - );  Run  L4  ( - ). 


By  contrast,  the  shear-layer  mode  shows  much  smaller  pressure  variations,  and 
mean  flow  streamlines  are  nearly  horizontal  along  the  mouth  of  the  cavity.  The 
low  pressures  (Cp  k  —0.08)  correspond  to  the  center  of  the  recirculation  region 
that  exists  in  the  rear  two-thirds  of  the  cavity  (see  figure  3.4).  The  low  pressure 
region  appears  to  be  the  result  of  the  vortical  swirling  motion  in  the  cavity,  and 
not  because  of  an  expansion  of  the  flow  into  the  cavity  as  was  the  case  for  wake 
mode.  There  is  a  very  small  impingment  region  at  the  rear  edge,  where  Cp  reaches 
about  0.1.  However,  over  the  majority  of  the  rear  face  Cp  <  0,  which  is  consistent 
with  the  measurements  of  Roshko  [1955]  and  others. 

Similar  flow  features  are  also  evident  in  the  mean  streamwise  velocity  profiles, 
shown  in  figure  3.18.  Wake  mode  oscillations  have  a  substantial  effect  on  the  mean 
even  as  far  as  1-2  depths  above  the  cavity,  and  also  affect  the  upstream  boundary 
layer,  to  about  a  quarter  of  a  depth  upstream.  As  we  noted  in  section  3.1.2,  the 
nominal  momentum  thickness  that  we  use  is  that  of  the  initial  condition,  and 
figure  3.18  shows  that  this  is  indeed  substantially  modified  by  the  flow  in  wake 
mode.  It  is  modified  only  by  a  very  small  amount  in  shear-layer  mode. 

In  Gharib  L  Roshko  [1987],  transitions  between  non-oscillatory,  shear-layer, 
and  wake  modes  occurred  at  L/9o  =  80  and  L/9o  =  160,  respectively.  The  present 
data  indicate  a  change  from  wake  mode  to  shear-layer  mode  that  depends  also  on 
the  Mach  number  and  Reynolds  number.  The  specific  parametric  dependence  of 
the  transition  to  w'ake  mode  is  discussed  in  section  3.3.3.  The  drag  on  the  cavity 
(given  for  selected  runs  in  table  3.1),  is  of  comparable  magnitude  to  the  values 
reported  by  Gharib  &  Roshko  [1987]  for  the  different  regimes.  An  impingment  of 
the  flow  on  the  rear  step  was  also  noted  in  the  experiments,  and  it  also  appears  (see 
Figure  7(a)  of  Gharib  L  Roshko  [1987])  from  dye  visualizations  that  the  boundary 
layer  separates  upstream  of  the  cavity  leading  edge. 
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3.3.2  Frequency  of  vortex  shedding 

The  spectra  of  the  resonant  instabilities  in  wake  mode  are  very  different  from 
those  of  the  shear- layer  mode  (section  3.2.2).  After  an  initial  transient,  which 
at  early  times  is  similar  to  shear- layer  mode,  the  flow  becomes  nearly  periodic 
in  time,  with  the  fundamental  period  corresponding  to  the  vortex  shedding  from 
the  leading  edge  (see  figure  3.16).  The  spectrum  (after  the  onset  of  wake  mode) 
consists  of  a  dominant  frequency  and  strong  peaks  at  its  harmonics.  The  4M  series 
of  runs  transitioned  from  shear-layer  to  wake  mode  between  0.3  <  M  <  0.4,  and 
the  peak  frequencies  were  plotted  in  figure  3.9.  For  M  >  0.3,  the  peaks  fall  well 
below  the  first  Rossiter  mode  prediction,  and,  unlike  the  shear-layer  mode,  they 
show  little  variation  with  M.  For  M  =  0.4  to  M  =  0.8  the  fundamental  frequency 
varied  less  than  4%,  compared  to  the  expected  variation  of  about  20%  for  Rossiter 
Mode  I.  The  4%  variation  is,  in  fact,  within  the  uncertainty  associated  with  the 
total  sampling  period  used  to  determine  the  frequency. 

The  lack  of  variation  with  M  indicates  that  w^ake  mode  is  not  a  result  of 
acoustic  feedback,  and  it  appears  that  the  feedback  in  this  case  is  provided  by 
the  complicated  recirculating  flow  in  the  cavity.  This  is  discussed  further  in  sec¬ 
tion  3.3.4. 

3.3.3  Parametric  dependence  of  the  transition  to  wake  mode 

In  table  3.1.  we  have  indicated  the  dominant  mode  (steady  flow,  shear-layer  mode, 
or  wake  mode)  for  the  different  runs,  and.  for  the  wake  mode  cases,  the  shedding 
frequency,  scaled  with  the  cavity  depth  and  freestream  velocity.  Once  established, 
the  shedding  frequencies  in  wake  mode  are  all  nearly  the  same  for  cavities  ■with 
L/D  =  4,  meaning  that,  in  addition  to  the  invariance  to  Mach  number  noted 
above,  they  are  not  influenced  by  L/6q  or  D/6q.  A  lack  of  dependence  of  the 
shedding  frequency  on  the  boundary-layer  thickness  is  a  common  feature  of  vortex 
shedding  behind  bluff  bodies  (for  instance,  for  laminar  flow  over  a  cylinder,  it  is 
well  known  that  the  shedding  frequency  St  —  fD/U  is  constant  over  a  wide  range 
of  Reynolds  numbers).  We  note  that  for  the  cavity  with  LjD  —  5  (run  L5),  the 
shedding  frequency  is  about  10%  lower  than  that  of  the  L/D  =  A  runs,  and  that  a 
periodic  wake  mode  was  not  detected  in  any  of  the  runs  with  L/D  —  3  or  smaller. 
Several  of  the  L/D  =  3  runs  exhibited  some  characteristics  of  both  shear-layer 
mode  and  wake  mode. 

While  the  wake  mode  tj'pically  occurs  for  larger  L/6q  and  D/9q,  the  transition 
depends  as  well  on  the  Mach  and  Rejmolds  numbers.  For  example  the  runs  L4 
and  TK4b  both  have  L/D  —  A,  but  the  latter  is  oscillating  in  shear-layer  mode,  as 
shown  by  the  instantaneous  vorticity  fields  plotted  in  figure  3.19,  which  should  be 
contrasted  to  figure  3.15.  Moreover,  the  transition  to  w^ake  mode  also  depends  on 
M ,  but  it  is  interesting  that,  as  previously  noted,  once  wake  mode  is  established 
the  frequency  of  vortex  shedding  is  nearly  independent  of  M.  The  4M  series  of 
runs  shows  transition  to  wake  mode  as  the  Mach  number  is  increased,  holding 
other  parameters  constant.  The  TK4  series  of  runs  shows  evidence  of  transition  to 
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Figure  3.19:  Instantaneous  vorticity  contours  for  run  TK4b  (shear-layer  mode). 
Contour  levels  same  as  figure  3.15. 


wake  mode  as  the  D/6q  is  increased,  holding  L/D,  M,  and  Re$  nearly  constant, 
and,  finally,  the  R4  series  shows  transition  with  increasing  Ree  when  L/D,  M.  and 
D/6q  are  held  nearly  constant. 

In  all  of  the  simulations  that  have  been  performed,  examination  of  the  instanta¬ 
neous  flow  fields  shows  that  the  initial  growth  of  instability  is  always  in  shear-layer 
mode,  where  acoustic  radiation  from  the  downstream  edge  is  the  mechanism  for 
the  feedback.  The  shear-layer  oscillations,  at  early  times,  are  growing  in  amplitude 
with  each  subsequent  cycle,  and  the  wake  mode  sets  in  only  after  the  oscillations 
have  reached  a  certain  threshold.  Thus  it  appears  that  the  emergence  of  wake 
mode  is  determined  by  the  conditions  in  shear-layer  mode,  and  this  forms  the  ba¬ 
sis  of  our  attempt  below  to  model  the  parametric  dependence  of  the  wake  mode. 
\^e  note  that  we  have  not  attempted  to  look  for  hysteresis  in  the  transition  to 
wake  mode,  as  there  is  no  unambiguous  way  to  do  so  computationally. 

3.3.4  Convective  and  absolute  instability 

As  indicated  above,  the  oscillaticu  frequency  in  wake  mode  is  independent  of  the 
-Mach  number,  over  a  large  range  from  0.4  <  M  <  0.8.  This  implies  that,  unlike 
the  oscillations  described  by  the  Rossiter  mechanism,  wake  mode  oscillations  are 
not  acoustically  forced  convective  instabilities.  In  wake  mode,  the  flow  exhibits 
several  qualitative  features  of  absolutely  unstable  parallel  flows,  as  described  in, 
e.g..  Huerre  k  .Monkewitz  [1985],  and  this  may  provide  an  explanation  of  the 
governing  mechanism. 

The  notion  of  absolute  or  convective  instability  is  a  useful  tool  for  describing 
self-sustained  oscillations  in  parallel  shear  flows.  The  idea  is  the  following:  given  a 
certain  parallel  flow,  if  energ>^  from  a  disturbance  grows  in  time,  but  travels  only 
downstream,  the  flow  is  called  convectively  unstable.  In  such  flows,  there  is  no 
possibility  for  self-sustained  oscillations,  since  if  the  flow  disturbance  is  removed, 
the  oscillations  convect  away.  However,  if  energy  grows  and  travels  upstream, 
then  the  flow  is  called  absolutely  unstable,  and  such  flows  have  the  potential  for 
self-excitation.  In  this  context,  flows  that  exhibit  self-sustained  oscillations  are 
called  globally  unstable  [Huerre  k  Monkewitz,  1990].  This  terminology^  is  normally 
reserved  for  (nearly)  parallel  flows,  however— non-parallel  flows  will  normally  not 
be  steady  solutions  of  the  Navier-Stokes  equations,  so  it  does  not  make  sense  to 


3.3  Wake  mode 


43 


discuss  their  stability. 

For  cavity  oscillations  in  shear-layer  mode,  the  flow  is  nearly  parallel,  so  these 
stability  concepts  apply  directly.  The  established  view  of  cavity  oscillations  in 
shear-layer  mode  is  one  of  a  convective  Kelvin-Helmholtz  instability  in  the  shear 
layer,  which  leads  to  global  instabihty  through  pressure  feedback  from  acoustic 
waves  generated  near  the  trailing  edge.  Indeed,  if  D/Oq  ':$>  so  that  the  boundary 
condition  at  the  cavity  bottom  may  be  neglected,  and  if  flow  within  the  cavity 
is  relatively  quiescent  (such  that  the  velocity  ratio  across  the  shear  layer  is  ap¬ 
proximately  unity),  then  analysis  shows  that  the  shear  layer  is  indeed  (locally) 
convectively  unstable  [Huerre  Monkewitz,  1985].  The  good  prediction  of  oscilla¬ 
tion  frequencies  by  the  Rossiter  formula,  and  the  computational  confirmation  given 
earher  (e.g.,  the  results  of  section  3.2.4),  confirm  the  feedback  process  leading  to 
global  instability. 

In  wake  mode,  the  flow  past  the  cavity  is  far  from  parallel,  and  the  distur¬ 
bances  are  far  from  small,  so  the  notions  of  absolute  and  convective  stability  do 
not  strictly  apply.  Nevertheless,  these  notions  have  been  applied  to  very  similar 
oscillating  wake  flows,  such  as  the  wake  behind  a  cylinder,  and  the  predictions  of 
the  theory  have  been  shown  to  be  in  good  agreement  with  experiment,  despite 
these  non-parallel,  large-amplitude  eS*ects  [Huerre  &  Monkewitz,  1990].  We  there¬ 
fore  apply  these  ideas  to  the  cavity  flow  in  wake  mode,  even  though  the  theoretical 
justification  is  questionable. 

With  these  ideas  in  mind,  the  flow  in  wake  mode  is  also  globally  unstable,  since 
self-sustained  oscillations  persist  despite  no  external  forcing.  However,  since  the 
frequency  is  independent  of  Mach  number,  the  feedback  mechanism  is  apparently 
no  longer  acoustic.  It  is  therefore  plausible  that  absolute  instability  provides  the 
mechanism  necessary  for  global  instability. 

The  mean  flow  profiles  shown  in  figure  3.18  show  that  in  shear  layer  mode 
(run  L2).  there  is  very  little  backflow  in  the  cavity,  and  furthermore  the  backflow 
is  confined  to  a  region  at  the  rear  of  the  cavity.  By  contrast,  the  wake  mode 
profiles  (run  L4)  reveal  a  substantial  backflow  even  very  close  to  the  cavity  leading 
edge.  Huerre  k:  Monkewitz  [1985]  have  shown  that  for  tanh  profiles  with  greater 
than  13.69c  backflow,  the  flow  is  absolutely  unstable.  The  backflow  in  the  mean 
profiles  for  run  L4  reaches  a  maximum  of  over  389c  of  the  freestream  velocity 
near  the  center  of  the  cavity,  decreasing  to  about  21%  within  two  bound  ary- layer 
thicknesses  of  the  leading  edge,  and  15%  within  one  boundary-layer  thickness.  The 
backfiow  for  run  L2  is  much  smaller,  reaching  a  maximum  of  23%  of  the  freestream 
near  the  rear  of  the  cavity,  and  never  exceeding  3%  of  the  freestream  anyv^here  in 
the  front  40%  of  the  cavity.  The  mean  profiles  from  run  L4  are  not  well  described 
by  tanh  profiles,  but  the  main  point  is  that  there  is  significant  backflow  in  the 
cavity,  and  in  this  sense  the  profiles  are  qualitatively  similar  to  profiles  which  have 
been  shoum  to  be  absolutely  unstable.  It  is  possible,  then,  that  absolute  instability 
may  provide  a  mechanism  for  transition  to  wake  mode.  It  is  also  conceivable  that 
similar  ideas  to  those  expressed  in  Monkewitz  k  Nguyen  [1987]  might  be  used  to 
predict  the  shedding  frequency  in  wake  mode. 
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3.4  Prediction  of  the  oscillation  regimes 

Here  we  discuss  a  method  of  predicting  transition  between  the  various  flow  regimes: 
DO  oscillations,  shear-layer  mode,  and  wake  mode.  Our  goal  is  not  a  quantitative 
prediction  for  the  amphtude  of  oscillation,  but  rather  an  approximate  scaling  law 
to  determine  the  parametric  dependence  of  the  transitions  between  the  dififerent 
regimes. 

Our  criteria  for  predicting  transition  to  wake  mode  are  based  on  the  follow¬ 
ing  observations:  as  certain  parameters  (e.g.,  L/Oq)  are  increased,  the  Kelvin- 
Helmholtz  disturbances  in  the  shear  layer  grow  to  larger  amplitude.  The  larger 
amplitude  disturbances  induce  a  larger  recirculating  flow  within  the  cavity,  pos¬ 
sibly  generating  a  region  of  absolute  instabihty  in  the  shear  layer,  and  ultimately 
inducing  large-scale  vortex  shedding  from  the  leading  edge.  The  key  ideas  here  are 
that  larger  amplitude  disturbances  lead  to  larger  back  flow  in  the  cavity,  which 
ultimately  leads  to  wake  mode. 

The  onset  of  fluctuations  (i.e.,  the  transition  from  steady  flow  to  shear-layer 
mode)  demonstrates  a  similar  dependence  on  the  parameters  of  the  problem, 
since  the  amphfication  by  the  shear  layer  must  exceed  a  certain  threshold  for 
self-sustained  oscillations  to  occur  [Woolley  k  Karamcheti,  1974].  At  low  Mach 
number.  Sarohia  [1975]  determined  that  there  is  a  minimum  cavity  length,  relative 
to  the  (laminar)  incoming  boundary-layer  thickness,  for  which  cavity  oscillations 
may  occur.  For  sufficiently  large  D/Oq  {D/Oq  >~  15)  he  found  oscillations  only 
when 


\/Ree  L/6o  >~  800.  (3.4) 

which  also  shows  that  there  is  a  minimum  speed,  a  maximum  viscositv,  and  a 
maximum  boundary-layer  thickness,  holding  other  parameters  constant,  beyond 
which  there  are  no  oscillations.  Gharib  k  Roshko  [1987]  yJso  measured  the  onset 
of  fluctuations  at  ^/We  L/Oq  %  780.  As  D/Oo  was  decreased  below  about  15, 
Sarohia  [1975.  found  a  rapid  increase  in  the  minimum  length  for  oscillations  to 
occur.  The  data  for  turbulent  fluctuations  and  higher  \Iach  numbers  are  not  so 
comi)lete.  but  the  general  trend  is  that  the  minimum  length  increases  with  tur¬ 
bulence  [Krishnamurty,  1956:  Sarohia.  1975],  and  decreases  with  increasing  Mach 
immber  [Krishnamurty,  1956]. 

Thus,  if  the  overall  amplification  of  disturbances  can  be  predicted,  simple  cri¬ 
teria  for  both  the  onset  of  fluctuations  and  the  transition  to  wake  mode  would  be 
that  the  amplification  exceeds  a  certain  value.  In  order  to  test  this  hypothesis,  we 
estimate  the  amplification  based  on  ideas  about  the  shear-layer  instability,  and  the 
efficiency  of  the  acoustic  radiation  from  the  trailing  edge.  Although  some  of  the 
^PP^^^imations  made  are  crude,  our  goal  is  to  establish  a  rough  parametric  depen¬ 
dence  of  transitions  between  the  flow  regimes,  rather  than  a  detailed  prediction 
for  the  overall  amplitude. 
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Figure  3.20:  Simple  model  of  nonlinearities. 

3.4.1  The  role  of  nonlinearities 

All  our  analysis  of  the  governing  physical  mechanisms  will  be  linear.  One  might 
object  at  the  outset  that  linear  amplification  rates  say  nothing  about  the  final 
amplitude  of  a  limit  cycle — in  fact,  once  the  fiow  has  settled  into  a  limit  cycle, 
the  loop  gain  (total  amplification  once  around  the  feedback  loop)  must  be  unity, 
regardless  of  the  amplitude  of  the  limit  cycle  (otherwise,  oscillations  would  grow  or 
decay).  One  must,  in  fact,  take  nonlinearities  into  account  to  predict  amplitudes. 

However,  under  simple  yet  reasonable  assumptions  about  the  nonlinearities 
present,  the  linear  part  of  the  loop  gain  defined  above  does  predict  both  the  onset 
of  oscillations  and  the  final  amplitude  of  a  limit  cycle.  Assume  the  system  may  be 
approximated  (crudely)  by  linear  dynamics  coupled  with  a  frequency-independent 
nonlinear  saturation,  as  shown  in  figure  3.20.  Here,  G{s)  is  the  transfer  function 
of  the  linear  elements  in  the  model  (for  our  purposes,  this  will  be  the  shear-layer 
amplification  and  scattering  at  the  downstream  comer),  and  %!;{')  is  a  saturation 
nonlinearity,  defined  precisely  below.  This  form  of  the  model  lends  itself  to  de¬ 
scribing  function  analysis  [e.g.,  Khalil,  1996],  a  simple  tool  for  approximating  both 
frequencies  and  amplitudes  of  limit  cycles  in  nonlinear  systems. 

Describing  functions 

If  ?/  1.*=^  periodic,  of  the  form  y  =  AsinWt,  then  t^(y(t))  will  also  be  periodic,  with 
the  same  frequencN'.  and  may  be  wntten  as  a  Fourier  series 

oc 

=  y^^Ck{A)smkijt, 
k=i 

where  we  have  assumed  that  if  is  odd.  In  the  simplest  form  of  describing  function 
analysis,  one  neglects  higher  harmonics  (assuming,  for  instance,  that  they  will  be 
attenuated  by  G(s)),  and  considers  only  k  =  1.  defining  the  describing  function 

N(A)  —  [  if{Asmt)smtdt. 

A  ttAJo 


(3.5) 
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Sinusoidal  solutions  y  —  ylsinwt  of  the  feedback  loop  must  then  satisfy  y  — 
G{i(jj)'ip{y)  ss  G{i(ju')N{A)y,  which  gives  the  harmonic  balance  equation 

G{iuj)N{A)  =  1.  (3.6) 

For  odd  nonlinearities,  N{A)  is  real,  and  so  taking  the  phase  of  this  equation  yields 

ZG(iw)  =  0,  (3.7)  '■ 

which  determines  possible  frequencies  of  oscillation.  (This  is  the  same  criterion, 
used  in  Rossiter’s  model  for  the  resonant  frequencies.)  Taking  the  magnitude 
of  (3.6)  gives 


|G(zu;)|  =  l/NiA),  (3.8) 

which  determines  the  amplitude  A  of  the  limit  cycle  with  frequency  u.  Stability  of 
the  limit  cycle  may  be  addressed  separately,  using,  for  instance,  Nyquist  stability 
techniques. 

The  key  assumption  in  the  above  discussion  is  that  one  may  neglect  higher 
harmonics  {k  >  1).  For  the  cavity  flow,  the  Kelvin-Helmholtz  mechanism  in  the 
shear  layer  amplifies  low  frequencies,  but  attenuates  higher  frequencies,  dropping 
off  quite  sharply  as  the  frequency  increases  (see  figure  6.7).  Thus,  the  above 
assumptions  are  reasonable  in  the  present  context.  Note  that  rigorous  results 
are  available,  and  one  may  rigorously  prove  existence  (or  non-existence)  of  limit 
cycles,  if  the  linear  part  satisfies  certain  conditions  [Khalil,  1996],  These  conditions 
quantify  how  sharply  the  linear  part  must  drop  off  as  frequency  increases.  Here, 
since  we  do  not  know  the  precise  nature  of  the  nonlinearities,  or  quantitative 
information  about  the  linear  dynamics,  we  cannot  directly  verify  these  conditions, 
but  they  are  not  difficult  to  check,  given  a  quantitative  model. 

Saturations 

We  now  prove  some  properties  of  describing  functions  for  saturation  nonlinearities, 
defined  as  follows: 

Definition  3.1  An  odd  function  v  is  a  saturation  nonlinearity  if  p{y) /y  is  positive 
and  decreasing  for  all  y  >  0,  with  'C^{y)/y  — ^  0  as  y  oc. 

Note  that  v(y)/y  need  not  be  strictly  decreasing.  We  now  demonstrate  some 
properties  of  describing  functions  for  this  class  of  nonlinearities. 

Proposition  3.1  Letip  be  a  saturation  nonlinearity.  Then  its  describing  function 
N  satisfies  the  following  properties: 

1.  A"(i4)  >  0  for  all  >1  >  0. 

2.  A'(i4)  0  as  i4  oc. 
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3.  N{A)  i/;'(0)  as  A  0. 

4-  N{A)  is  decreasing  for  A  >  0. 

Proof:  Property  1  follows  directly  from  equation  (3.5),  since  ilfy)  >  0  for  y  >  0, 
thus  making  the  integrand  positive  for  all  t.  Property  2  al&o  follows  from  (3.5), 
since  ii){y)/y —*  0  as  y oo\ 

2  t^(Asint)  9 

lim  N  (A)  =  lim  -  /  ^  sin^  t  dt 

/l— oo  >4^00  n  Jq  a  sm  t 


=  -  r  h 


ip{Asmt)  .  2.  j. 

lim  ,  - sin  tdt 

1—00  A  sm  t 


The  proof  of  property  3  is  similar.  Letting  K  =  ip'iO)  =  lim/,_o  t^(h)/h,  we  have 

lin,  mA)  ^  ^ 

A^o  A^orcJo  A  Sint 

2 

—  —  I  K  sin^  t  dt 


Property  4  uses  that  'ip{y)/y  is  decreasing  for  y  >  0.  Let  Ai,A2  G  M  with  0  < 
-4]  <  An-  Then  for  t  G  (0,7r), 

A2smt  >  i4i  sint  >  0 

7p{A2S\nt)  ^  ^^{Aismt) 

^2sint  ^  ^isint 

tl>(A2S\nt)  .  ^  ^  ipiAisint)  .  ^ 

=>  --i - i  sm  t  <  — — : - ^  sin  t 

A2  Ai 


■  [  ^iA2  si 
:  Jo 


sin  t)  sin  tdt  < 


> 

sin  t)  sin  t  dt 

n  Jo 


N{A2)  <  N{Ail 


and  hence  -V{.4)  is  decreasing  for  ^  >  0. 


These  properties  show  that  for  saturation  nonlinearities,  the  describing  function  A^(^) 
must  hav^e  the  general  form  depicted  in  figure  3.21, 

The  figure  shows  the  amplitude  A  on  the  vertical  axis,  and  the  reciprocal  of 
the  describing  function  on  the  horizontal  axis.  Recall  that  a  limit  cycle  with 
amplitude  A  must  satisfy  (3.8),  so  when  the  linear  gain  |G(2cj)|  =  1/N{A),  the 
amplitude  of  the  limit  cycle  (as  predicted  by  describing  function  analysis)  is  A,  as 
indicated  in  the  figure.  The  properties  of  saturation  nonlinearities  guarantee  the 
qualitative  features  of  the  curve:  specifically,  it  must  always  “bend  to  the  right,” 
since  1/A^(i4)  is  an  increasing  function  of  A.  and  obviously  the  amplitude  A  can 
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Figure  3.21:  Describing  function  N{A)  for  a  saturation  nonlinearity  ^(•).  Here 
K  —  the  gain  of  ip{-)  for  small  amplitudes. 

never  decrease  if  |G(iw)|  increases.  (Otherwise,  N{A)  would  have  to  be  multiple¬ 
valued.)  Stability  of  the  limit  cycle  must  be  addressed  separately,  but  the  standard 
argument  based  on  the  Nyquist  criterion  (see  Khalil  [1996])  shows  that,  as  long 
as  the  approximations  in  describing  function  analysis  are  valid,  the  limit  cycle  is 
stable. 

Application  to  cavity  oscillations 

Without  knowing  details  of  the  nonlinearity  and  exact  vailues  for  G{iLj),  we  cannot 
predict  the  exact  value  of  the  amplitude.  But  if  we  obtain  scaling  laws  for 
the  amplitude  will  follow  (at  least  qualitatively)  the  same  scaling  laws:  for  |G(itj)! 
below  a  certain  threshold,  we  expect  no  oscillations  to  occur,  and  for  iG(ft.j)| 
above  this  threshold,  we  expect  a  stable  limit  cycle  whose  amplitude  increases 
with  |G(2lj)|. 

The  precise  way  in  which  nonlinearities  enter  into  the  physics  is  still  not  clearly 
understood.  One  scenario,  most  likely  to  occur  for  short  cavities  (small  Lf6o), 
is  that  employed  bj’  Cain  et  al.  [1996],  where  the  satmation  occurs  because  of 
spreading  of  the  shear  layer:  the  spreading  of  the  mean  flow  is  caused  by  Reynolds 
stresses  that  are  proportional  to  the  square  of  the  amplitude  of  the  oscillations,  plus 
any  background  turbulence  or  other  nonlinear  interactions  between  the  modes,  as 
described  by  Morris  et  al.  [1990],  Linear  growth  rates  of  the  Kelvin-Helmholtz  in¬ 
stability  waves  decrease  as  the  shear  layer  spreads,  and  eventually  become  negative 
for  very  thick  shear  layers,  thus  providing  a  saturation  mechanism.  For  this  mech¬ 
anism.  then,  the  shear-layer  disturbances  are  well  described  by  linear  methods, 
given  the  correct  shear-layer  thickness.  Our  data  in  section  3.2.4  are  consistent 
with  this,  as  shown  by  the  good  agreement  of  the  data  with  linear  predictions 
(see  figure  3.12).  Though  the  amplification  predicted  by  hnear  theory  does  not 
agree  well  (figure  3.14),  this  is  more  likely  because  of  viscous  effects,  rather  than 
nonlinear  effects,  bec8,use  there  is  disagreement  even  for  small  x/0o. 

The  saturation  mechanism  more  commonly  considered  is  that  once  the  oscil¬ 
lations  grow  large,  the  linearization  is  no  longer  vahd,  and  the  resulting  nonlinear 
effects  limit  the  growth.  This  effect  is  most  likely  for  longer  shear  layers  (large 
L/Oo),  as  in  the  experiments  of  Cattafesta  et  al.  [1997],  which  had  a  very  thin 
upstream  boundary  layer  (L/0o  =  328).  Knisely  &  Rockwell  [1982]  also  attributed 
their  saturation  to  this  mechanism,  although  it  is  also  plausible  that  the  saturation 
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might  be  explained  by  mean  flow  spreading,  as  the  linear  stability  calculation  they 
compared  to  had  a  constant-thickness  shear  profile. 

Both  of  these  forms  of  nonlinearities  fit  nicely  into  the  describing  function 
framework  described  above;  if  a  shear-layer  disturbance  at  the  leading  edge  is 
given  by  f(0,t)  =  Aexp(iujt)  +  c.c.,  then  at  a  distance  x  downstream,  the  amplified 
disturbance  is  approximately  described  by 

f(x,t)  =  Nx(A)Aexp(iiJt)  ■  Gx(iuj)  +  c.c.,  (3.9) 

where  Gx{iuj)  is  the  linear  gain  and  phase  at  location  x  (for  instance,  Gx(iu;)  = 
exp(2a(a;)x)),  and  NxiA)  is  a  nonlinear  saturation,  equal  to  unity  for  small  A,  and 
decreasing  for  large  A,  and  also  depending  on  the  position  x.  Thus,  NxiA)  in  (3.9) 
plays  the  same  role  as  the  describing  function  above,  and  similar  conclusions  about 
the  stability  and  amplitude  of  the  limit  cycle  follow. 

3.4.2  Shear  layer 

The  shear-layer  amplification  is  a  complicated  function  of  M,  L/6o,  D/6q,  and 
Reg.  To  proceed,  we  perform  inviscid  linear  stability  calculations  similar  to  those 
described  in  section  3.2.4  for  each  set  of  computational  parameters  in  table  3.1,  and 
for  each  of  the  first  three  Rossiter  frequencies  (for  which  we  use  equation  (3.2)).  We 
assume  a  tanh  velocity  profile  with  a  linear  spreading  rate  dS^/dx  =  0.1.  Based 
on  the  discussion  above  and  in  section  3.2.3,  it  is  clear  that  the  true  spreading 
rate  is  dependent  on  growth  of  the  disturbances  in  the  layer,  which,  in  turn,  is 
dependent  on  all  of  L/do,  D/Oq,  M,  and  Reg.  However,  we  performed  calculations 
using  different  values  of  the  spreading  rate,  and  found  the  amplitude  variation  to 
be  relatively  weak  compared  to  the  other  factors. 

3.4.3  Scattering 

In  order  to  assess  the  efficiency  of  the  sound  generation  at  the  edge  as  a  function 
of  the  Mach  number,  we  use  the  simple  model  proposed  by  Tam  Block  [1978], 
where  the  process  is  idealized  as  an  oscillating  compact  mass  source  (monopole). 
We  determine  the  strength  Q  of  the  monopole  source  by  setting  Q  equal  to  the 
mass  flow  rate  out  of  the  cavity.  Tam  points  out  that  “the  unsteady  mass  addi¬ 
tion  and  removal  at  the  trailing  edge  of  the  cavity  as  the  cause  of  the  acoustic 
disturbance  . . .  essentially  suggests  a  dipole  source  at  the  trailing  edge,"  since  by 
this  mechanism  a  compression  wave  outside  the  cavity  (mass  addition)  would  cor¬ 
respond  to  a  rarefaction  wave  inside  the  cavity  (mass  removal).  However,  as  Tam 
points  out,  experimental,  and  now  computational,  observations  indicate  that  the 
pressure  waves  inside  and  outside  the  cavity  are  in  fact  in  phase,  which  corresponds 
to  a  monopole  source  at  the  trailing  edge.  Thus,  we  model  the  acoustic  source  as 
a  monopole  at  the  trailing  edge,  whose  strength  is  determined  by  the  mass  flow 
rate  out  of  the  cavity. 

If  density  fluctuations  are  small,  the  mass  flow  rate  per  unit  volume  is  given 
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by 


Q^rn  ^ 

^  V  V 


(3.10) 


where  v  is  the  vertical  velocity,  and  V  is  the  volume  of  the  source,  at  this  time  un¬ 
known,  Assuming  the  velocity  perturbations  are  sinusoidal,  with  an  exponentially 
growing  envelope,  the  integral  above  will  be  dominated  by  the  portion  closest  to 
the  trailing  edge,  and  thus  scales  as  A|i;(L)|,  where  \v{L)\  denotes  the  amplitude 
of  the  velocity  perturbation  near  the  traihng  edge,  and  A  is  the  wavelength  of  the 
instability  wave.  Expressing  the  wavelength  as  A  =  27TCp/u,  where  Cp  is  the  phase 
speed  of  the  instability  wave,  the  source  strength  scales  as 


Q  =  iujQ  oc  iuj 


PoCp\v{L)\  f)oCp\v{L)\ 


OC 


u 


7 


(3.11) 


This  expression  may  also  be  arrived  at  on  dimensional  grounds  alone,  if  it  is 
assumed  that  the  monopole  strength  is  proportional  to  the  vertical  velocity  (or 
shear-layer  displacement)  at  the  trailing  edge. 

Density  fluctuations  then  satisfy  a  Helmholtz  equation  [see,  e.g.,  Crighton, 
1975],  and  far  from  the  source  are  given  by 


P  Cp\v{L)\  e*''’’ 


(3.12) 


where  r  is  the  distance  from  the  source,  and  k  —  uj/a^c  is  the  wavenumber.  Note 
that  there  are  two  asymptotic  arguments  leading  to  equation  (3.12).  The  first  is 
that  the  retarded  time  is  negligible  over  the  source  region  (compact  source),  and 
the  second  is  that  we  are  at  large  distance,  on  the  scale  of  the  wavelength,  from  the 
source.  The  second  approximation  is  not  valid  at  very  low  Mach  numbers,  where 
the  cavity  length  is  small  compared  to  the  wavelength  of  the  acoustic  waves.  For 
the  present  range  of  il/,  the  asymptotic  properties  of  the  Green’s  function  for  the 
Helmholtz  equation  (i.e.,  the  zeroth  order  Hankel  function  of  the  first  kind)  can 
be  used  to  show  that  there  is  negligible  error  in  the  second  assumption. 

At  the  cavity  leading  edge,  r  =  L,  we  then  have  the  following  scaling  law 
(assuming  Cp/U  =  const,  as  in  the  Rossiter  model): 


Ml  oc  (3.13) 

Poc  U 

where  ip(0)|  denotes  the  magnitude  of  density  fluctuations  at  the  leading  edge, 

X  =  0. 


3.4.4  Loop  gain 

We  now  have  enough  information  to  determine  the  scaling  of  the  loop  gain  for 
the  cavity;  that  is,  the  amplification  of  a  disturbance  as  it  is  amplified  by  the 
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Figure  3.22:  Correlation  of  loop  gain  0  with  mode  of  cavity  oscillations.  All  runs 
from  table  3.1  with  56  <  Re^  <  70  are  shown.  As  in  table  3.1,  NO  corresponds 
to  no  oscillations,  SL  denotes  shear-layer  mode,  M  denotes  mixed  mode,  and  W, 
wake  mode. 

shear  layer  and  converted  into  a  monopole  acoustic  source.  (We  do  not  include 
receptivity  effects  in  this  crude  model.)  Denoting  the  shear-layer  amplification  by 
A{St.M,L/6o.D/Oo),  the  loop  gain  scales  as 

0n  oc  A(5t„,  M.  L/do,  D/eQ)M^I^St-^>^  (3.14) 

where  St-n  denotes  the  Strouhal  number  of  Rossiter  mode  n,  given  by  Rossiter’s 
formula  (3.2).  The  total  amplification  is  then  d  =  maXn/3„,  and  we  take  /?  as 
our  criterion  for  predicting  the  onset  of  oscillations,  and  eventually  wake  mode. 
Note,  however,  that  the  dependence  of  Re\molds  number  is  not  included  in  (3.14), 
as  our  shear-layer  calculations  were  inviscid.  Variation  with  Re\molds  number  is 
considered  separately  below.  Note  also  that  this  scaling  law  is  not  valid  for  small 
M,  as  the  asymptotic  approximation  (3.12)  no  longer  holds. 

Figure  3.22  shows  how  (3  correlates  with  the  mode  of  oscillation.  Because 
Reynolds  number  effects  are  not  included,  we  compare  all  runs  in  table  3.1  that 
fall  in  a  limited  Refolds  number  range,  56  <  Ree  <  70.  For  each  run,  we  compute 
the  amplification  0n  for  the  first  three  Rossiter  modes  (n  =  1,2.3),  and  /?  for  each 
run  is  then  the  maximum  of  the  three.  Figure  3.22  shows  0  vs.  the  observed  mode 
of  cavity  oscillation,  and  w^e  find  that  for  this  Reynolds  number  range,  0  <  4.8 
corresponds  to  no  oscillations,  A.%  <  0  <  15  corresponds  to  shear-layer  mode, 
15  <  /?  <  21  for  mixed  mode,  and  0  >  2\  for  pure  wake  mode. 

3.4.5  Parametric  dependence 

We  now  discuss  some  results  related  to  the  specific  parametric  dependence  of  the 
transitions.  In  hguxe  3.23,  we  have  plotted  the  total  amplification  from  the  shear 
layer  alone  (corresponding  to  An  in  equation  (3.14))  for  the  first  three  Rossiter 
modes,  as  a  function  of  L/6o  and  DfOo,  at  M  =  0,6.  Since,  for  fixed  Mach 
number,  the  shear  layer  amplification  is  the  predominant  component  of  the  loop 
gain  /3,  the  figure  also  indicates  trends  in  /?  as  the  cavity  length  is  varied.  The 
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Figure  3.23:  The  shear-layer  amplification,  An,  versus  L/Bq  for  different  D/Bq. 
with  M  =  0.6.  The  dominant  mode  is  shown  at  each  L/Bq,  with  n  =  1  (o),  n  =  2 
(□),  and  n  =  3  (A). 

amplification  shows  strong  dependence  on  D/Bq,  especially  for  D/Bq  <~  20.  This 
very  well  explains  the  trend  seen  (for  low  Mach  number)  by  Sarohia  [1975],  who 
found  the  minimum  cavity  length  for  oscillations  increased  dramatically  for  small 
D ! Be .  The  plot  also  shows  why  the  higher  Rossiter  modes  tend  to  dominate  for 
longer  cavities,  as  has  been  observed  in  experiments. 

It  is  interesting  to  note  that  the  predictions  show  that  multiple  modes  are  likely 
to  be  present  at  roughly  equal  amplitudes  for  cavities  of  certain  length,  while,  for 
other  lengths,  one  mode  is  clearly  dominant.  For  example,  for  LjBe  =  125  and 
DjBe  —  30,  modes  2  and  3  shoi.dd  be  of  roughly  equivalent  strength  (ignoring, 
nonlinear  interactions  and  frequency  dependence  of  scattering/receptivity). 

Figure  3.23  also  shows  similar  effects  for  the  transition  to  wake  mode  as  a 
function  of  DjBe.  but  given  that  the  amplitude  threshold  is  higher,  wake  mode 
should  not  occur  in  cavities  with  DJBq  <~  20. 

In  figure  3.24  the  total  amphfication  in  the  shear  layer  for  the  first  three  Rossiter 
mode.s  is  plotted  as  a  function  of  LjBe  and  A/,  for  DjBe  =  30.  The  amplification 
rates  are  highest  for  the  lowest  Mach  numbers,  which  is  consistent  with  many 
pre^•ious  analyses  of  the  Kelvin-Helmholtz  instability  in  compressible  shear  layers. 
This  decreasing  amplitude  must  be  balanced  by  the  increase  in  radiation  efficiency 
at  the  edge  as  M  is  increased.  As  noted  above,  it  is  observed  in  experiments 
[Krishnamurty,  1956]  that  the  minimum  cavity  length  for  the  onset  of  fluctuations 
decreases  with  A/ .  This  is  indeed  what  the  model  predicts,  as  evidenced  by  the 
good  agreement  with  our  data  in  figure  3.22. 

\N  e  note  that  the  computational  data  dso  shows  a  tendency  for  transition  to 
wake  mode  with  increasing  Reynolds  number.  The  R4  series  of  runs,  for  example, 
show  a  change  from  shear-layer  to  wake  mode  for  LjBe  »  75  as  Re^  varies  from 
40  to  80.  The  model  of  the  transition  process  discussed  above  is  based  on  an 
inviscid  stability  analysis,  and  is  therefore  unable  to  predict  any  Reynolds  number 
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L/00 

Figure  3.24:  The  shear-layer  amplification,  versus  L/6o  for  different  A/,  with 
D/6q  =  30.  The  dominant  mode  is  shown  at  each  L/6o,  with  n  =  1  (o),  n  —  2 
(□),  and  n  =  3  (a), 

effects.  However,  for  the  range  of  low  Re}molds  numbers,  substantial  reduction 
in  the  amplitude  gain  in  the  sheeir  layer  may  be  expected  [e.g.,  Michalke,  1984], 
which  is  consistent  with  figure  3.14,  and  with  the  observed  trend  in  transition  to 
wake  mode. 


3.5  Concluding  remarks 

We  have  used  numerical  simulations  to  explore  the  operating  regimes  of  the  lami¬ 
nar  flow  past  a  two-dimensional  rectangular  cavity.  For  short  cavities  (relative  to 
the  upstream. boundary-layer  thickness),  and  for  low  Mach  numbers  and  Reynolds 
numbers,  the  flow  is  steady.  As  these  parameters  are  increased,  the  flow  transi¬ 
tions  into  a  shear-layer  mode,  where  self-sustained  oscillations  occur.  This  is  the 
regime  usually  observed  in  experiments,  and  the  acou'stic  fields  predicted  by  the 
simulation  agree  well  with  schlieren  photographs  from  the  experiments  by  Krish- 
namurty  [1956],  The  growth  of  disturbances  in  the  shear  layer  is  predicted  well  by 
a  locally-parallel  linear  stability  calculation,  where  the  thickness  of  the  shear,  layer 
is  measured  from  the  simulations. 

For  longer  cavities,  and  larger  Mach  and  Reynolds  numbers,  the  flow  transitions 
into  a  wake  mode.  Similar  flow  features  have  been  observed  in  axisymmetric  cavity 
experiments  by  Gharib  Roshko  [1987],  and  in  experiments  in  a  pipe  with  closed 
side  branches  [Kriesels  et  ah,  1995].  The  frequency  of  oscillations  in  wake  mode 
is  independent  of  Mach  number,  indicating  a  purely  hydrodynamic  (non-acoustic) 
instability.  In  wake  mode,  a  significant  backflow  is  present  inside  the  cavity,  and  we 
hypothesize  that  this  backflow  leads  to  an  absolute  instability,  which  may  provide 
the  feedback  mechanism  leading  to  wake  mode.  We  have  used  a  simple  linear 
model  to  predict  the  scaling  laws  governing  transition  between  the  flow  regimes, 
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and  the  model  agrees  well  with  data  from  our  simulations. 

To  our  knowledge,  wake  mode  oscillations  have  not  been  observed  in  experi¬ 
ments  in  two-dimensional  cavity  geometries,  for  a  similar  range  of  L/D  as  we  have 
investigated,  (However,  very  recent  three-dimensional  PIV  measurements  in  an 
L/D  =  5  cavity  show  a  time- averaged  flow  with  a  significant  recirculating  region, 
reaching  about  15%  of  the  freestream  velocity  [Krothapalli  2001].)  According  to 
the  criteria  proposed  in  the  previous  section,  transition  to  wake  mode  occurs  when 
the  amplitude  of  cavity  oscillations  increases  beyond  some  threshhold.  Thus,  any 
effect  that  decreases  this  amplitude  should  inhibit  transition  to  wake  mode.  For 
instance,  turbulent  upstream  boundary  layers  will  likely  inhibit  transition  to  wake 
mode,  since  Krishnamurty  [1956]  observed  that  the  acoustic  radiation  is  less  intense 
when  the  boundary  layer  is  turbulent.  Furthermore,  three-dimensional  effects  mav 
significantly  change  the  character  of  the  instabilities  in  the  shear  layer,  making 
them  less  coherent,  and  ultimately  inhibiting  transition  to  wake  mode.  It  appears 
that  these  three-dimensional  effects  do  not  affect  the  fundamental  mechanisms  un- 
derl\  ing  shear-layer  mode,  as  our  two-dimensional  simulations  in  shear-layer  mode 
agree  well  with  experiments,  but  three-dimensional  effects  may  reduce  the  overall 
amplitude  of  oscillations,  which  would  inhibit  transition  to  w^ake  mode. 

In  the  remainder  of  this  thesis,  we  focus  on  the  modeling  and  control  of  shear- 
layer  mode,  both  because  it  is  the  regime  usually  observed  in  experiments,  and 
because  the  underlying  mechanisms  are  better  understood,  and  lend  themselves 
more  readily  to  the  modehng  approaches  w^e  discuss  in  the  following  chapters. 


Chapter  4 

Model  Reduction  for  Fluids 


In  this  chapter,  we  discuss  basic  techniques  for  obtaining  reduced-order  models 
from  simulations  of  fluid  flows.  First,  we  give  an  overview  of  the  basic  tools  we 
use.  Proper  Orthogonal  Decomposition  and  Galerkin  projection.  We  briefly  outline 
how  these  tools  have  been  apphed  to  incompressible  flows  in  section  4.3.  The  main 
result  of  this  chapter  is  the  apphcation  to  compressible  flows  in  section  4.4. 

4.1  Proper  Orthogonal  Decomposition 

The  central  idea  of  the  Proper  Orthogonal  Decomposition  (POD)  is,  given  a  col¬ 
lection  of  data,  elements  of  a  vector  space,  to  find  a  subspace  of  fixed  dimension, 
which  is  ‘‘optimal,”  in  the  sense  that  the  error  in  the  projection  onto  the  sub¬ 
space  is  minimized.  Thus,  the  method  is  inherently  data-based,  and  provides  a 
way  of  distilling  the  information  contained  in  large  datasets,  as  are  produced  by 
simulations.  The  solution  of  the  optimization  problem  is  surprisingly  simple — it 
reduces  to  a  standard  eigenvalue  problem,  and  provides  an  orthonormal  basis  for 
the  optimal  sabspace. 

Here,  we  describe  the  decomposition  in  the  context  of  an  abstract  Hilbert 
space.  We  use  the  abstract  setting  both  for  clarity,  and  more  importantly,  because 
in  subsequent  sections  we  will  want  to  use  different  inner  products,  to  obtain 
different  notions  of  “optimality.” 

The  usual  presentation  of  the  method  is  in  the  stochastic  setting  (see  Loeve 
[1978,  sec.  37.5]),  where  the  “data”  is  regarded  as  a  random  variable.  The  presen¬ 
tation  here  treats  the  data  as  a  finite  ensemble,  as  it  always  is  given  in  practice. 
This  treatment  removes  many  of  the  technical  difficulties  of  the  theory,  since  one 
need  not  worry  about  subtleties  of  infinite-dimensional  operators,  yet  it  retains  all 
of  the  essential  features  relevant  to  model  reduction. 

4.1.1  The  main  theorem 

First,  we  state  the  goal  more  precisely.  Let  ^  be  a  Hilbert  space,  with  inner 
product  (•,•).  The  goal  is,  given  an  ensemble  of  data  G  H  \  k  —  I, .. .  ,m}, 
find  a  subspace  S  of  fixed  dimension  n  <  such  that  the  error  E{\\u  —  F!9u||)  is 
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minimized.  Here,  ||  •  (|  is  the  induced  norm  on  H,  Ps  is  the  orthogonal  projection 
onto  the  subspace,  and  E{-)  denotes  an  average  over  k.  The  data  {u''}  could  be 
thought  of  as  an  ensemble  of  many  different  experiments,  or  as  a  time  average, 
with  representing  snapshots  of  a  fvmction  u{t)  at  different  times  t  =  tk-  We 
use  the  notation  u  to  represent  the  entire  ensemble  {u*}.  The  average  £'(•)  may 
be  a  weighted  average,  and  all  we  will  assume  is  that  it  is  linear,  and  thus  may  be 
interchanged  with  the  inner  product. 

For  our  purposes,  the  space  H  will  consist  of  functions  on  some  spatial  do¬ 
main  fl  in  which  a  fluid  evolves,  e.g.,  H  =  L^(Q).  These  functions  may  be  vector¬ 
valued,  and  in  applications  we  will  explicitly  state  which  inner  product  we  use. 

Let  E  H\j  =  \,. . .  ,n}  be  an  orthonormal  basis  for  the  subspace  S.  Then 
the  orthogonal  projection  is  given  by 


Psu  =  '^{u,<fj)<fj  (4,1) 

j=i 

and  it  is  simple  to  show  that  minimizing  E{\\u  —  -PswH)  over  S  is  equivalent  to 
maximizing  £(^5^11),  or  maximizing  E{\  {u.ipj)  j^)  over  orthonormal  func¬ 
tions  The  desired  ^pj  are  thus  extremals  of  the  functional 

J[v?]  =  £:(|(u.v?>|2)-A(||<,?||2-1).  (4.2) 

At  tiiLs  point,  it  is  not  obvious  that  the  extremals  of  J  will  be  orthogonal,  but  we 
will  see  shortly  that  they  are. 

Before  finding  extremals  of  (4.2).  we  introduce  some  notation.  The  space  of 
all  continuous  linear  complex-valued  functionals  on  H  is  called  the  dual  space  of 
H .  and  is  denoted  H'.  For  any  vector  v  e  H,  its  dual  vector  v*  €  H*  is  given  by 
!••(•)  =  (-.v).  a  linear  functional  on  H.  Clearly,  v’  depends  on  the  inner  product. 
For  re//  and  a  €  H',  the  tensor  product  r®  a  ://—»//  is  the  linear  map  given 
by  {i'  ®  a  )i-  =  va{v)  for  any  V’  €  H.  In  particular!  for  v.w  e  H,  we  have 

{v  ®  =  uu'*(t:’)  =  V  (t'.  ui) . 

which  of  course  depends  on  the  inner  product.  (See  Abraham  et  al.  [1988]  for  the 
general  theory  of  tensors  on  linear  spaces.) 

With  this  notation,  we  may  rewrite  the  quantity  |  (u,  ip)  |-  as 

I  (^!  V’)  P  =  <p)  (‘r- u)  =  ((y^,  u)  u,  ’p)  =  {{u  ®  u*)p,  p) , 

and  so.  introducing  the  linear  operator  R  :  H  H  given  by 

R-Eiu^u").  (4.3) 

we  may  rewrite  (4.2)  as 


J[p]^{Rp.p)-X{\\pf-l). 


(4.4) 
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It  follows  easily  from  the  definition  that  R  is  self-adjoint  (that  is,  (R(p.  i^)  =  Ri)) 
for  all  £  H).  Now,  a  necessary  condition  for  ^  to  be  an  extremal  of  J  is  that 
the  first  variation  vanish: 

J[p  +  5rP\  =  0,  for  all  'tp  E  H. 

(5=0 


A 

dS 


We  have 


<5=0 


J[(p  -h  ^  {R{ip  +  d^ip),  (p  +  S'tp)  -  A(  ((^  +  StP,  (f  +  S-ip)  -  l) 

do  I 

=  {Rxp,  if)  +  {Rp,  Ip)  -  A((t/;,  p)  +  {p,  i^)) 

=  2  Re  [  {Rp  —  Xp,  xp)  ] 


<5=0 


where  in  the  last  equation  we  have  used  that  R  is  self-adjoint.  If  p  is  an  extremal, 
this  quantity  must  vanish  for  all  vEuriations  xp,  and  so  we  have 

Rp  =  Xp,  (4.5) 


or  p  must  be  an  eigenfunction  of  R,  Now,  we  see  from  (4.5)  by  taking  an  inner 
product  with  p  that 

X  =  E{\{u,p)\^^).  (4.6) 

Thus,  R  is  positive  semi-definite  (A  >  0).  and  the  functions  pj  which  maximize 
our  original  functional  pj)  \^)  are  the  eigenfunctions  corresponding  to 

the  largest  n  eigenvalues  of  R. 

We  now  examine  properties  of  the  eigenvalues  and  eigenfunctions  of  R.  It  is 
obvious  from  the  definition  of  R  that  its  range  is  finite-dimensional  (it  is  spanned  by 
the  {u^}).  Furthermore,  R  is  bounded  as  long  as  the  are  bounded,  and  hence 
R  is  compact  (or  completely  continuous).  For  compact,  self-adjoint  operators, 
spectral  theory  guarantees  at  most  a  countable  infinity  of  nonzero  eigenvalues 
[Ricsz  k  Sz.-Nagy,  1990],  and  since  the  range  of  R  is  finite-dimensional,  the  number 
of  nonzero  eigenvalues  is  actually  finite.  Furthermore,  since  R  is  self-adjoint,  the 
eigenfunctions  may  be  chosen  to  be  orthonormal.  We  also  mention  that  the  solution 
to  the  optimization  problem 

n 

max ^  {Rpj.pj)  for  orthonormal  pj 
j=i 

may  be  obtained  directly  from  functional  analysis,  without  recourse  to  the  varia¬ 
tional  argument  given  above  [Riesz  k  Sz.-Nagy,  1990,  sec.  93]. 

We  summarize  these  results  in  the  following  theorem. 

Theorem  4.1  (Proper  Orthogonal  Decomposition)  Given  a  linearly  inde¬ 
pendent  set  {u^  e  H  \  k  ,m},  the  subspace  S  of  dimension  n  <  m  xjuhich 

minimizes  £^(||u  —  Ps^ll)  kas  an  orthonormal  basis  {pj  G  |  j  =  1, . . .  .n}  vjhere 
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(fj  are  solutions  of 


Ripj  —  (4.7) 

where  R  =  E{u  ®  iz  )  and  Ai  >  A2  >  ■  •  •  >  >  0  are  the  largest  n  eigenvalues 

of  R.  The  eigenfunctions  ^pj  are  called  POD  modes. 

We  examine  some  properties  of  the  POD  in  section  4.1.2,  but  first  we  illustrate 
with  some  examples. 

Finite  dimensions.  Let  H  =  with  the  standard  inner  product  (x,  y)  —  Xiyi. 
where  we  use  Einstein’s  convention  of  summing  over  repeated  indices.  Suppose  the 
'‘data  is  given  as  a  collection  of  m  finearly  independent  vectors  eW\k  = 
L  . . .  ,  777},  and  the  average  £*(•)  is  just  an  arithmetic  mean 

1  ^ 

(4.8) 

A:=l 

For  x.y  €  Rp,  we  have  (x®y*)y  =  nyj,  the  usual  dyadic  product,  so  the  operator 
/?  :  RP  — >  RP  is  the  matrix  given  by 


1  ^ 

Rij  -=  £(u  O  u")ij  =  ^  V 

k=\ 

a  pxp  real  symmetric  matrix.  Equation  (4.7)  is  then  a  standard  eigenvalue  problem 
on  RP. 


Scalar-valued  functions.  Let  H  =  1],  the  set  of  complex-valued  square- 

integrable  functions  on  the  interval  [0, 1],  with  the  inner  product 

{u,v)  =  f  v{x)u{x)dx. 

Jo 

Then  for  u.v.  e  L?[Q,  l], 

{uZv  )^{x)  ■=  u{x)  {^.v)  =  f  u{x)v{y)i^{y)  dy. 

Jo 

so  the  eigeinalue  problem  (4.7)  becomes 

/  r{u{x)u{y))^{y)dy  =  \>p{x), 

a  Fredholm  integral  equation,  with  kernel  K{x,y)  =  E{u{x)u{y)). 

Vector-valued  functions.  Now  let  H  =  C(fi,  V),  the  space  of  continuous  func¬ 
tions  from  some  (spatial)  domain  f2  c  to  a  vector  space  V  =  C^.  Define  an 
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inner  product  on  H  by 


(u,v)  =  /  v^(x)i4u(x)  dx, 

Jq 

where  ^4  €  is  a  positive-definite  Hermitian  matrix.  Then 

(u  ®  V*)v>(x)  =  u(x)  {(f,  v)  =  /  U(x)v^(y)>l(^(y)  dy 

Jn 

so  the  eigenvalue  problem  (4.7)  becomes 

f  E{u{x)u'^{y))A<p{y)dy  =  Xtp{x).  (4.9) 

Jn 

4.1.2  Properties  of  the  POD 

In  this  section  we  discuss  several  important  properties  of  the  POD. 

Energy  captured 

There  is  a  physical  significance  to  the  value  of  the  eigenvalue  in  (4.7):  Ajt 
represents  the  average  “energy'”  captured  by  POD  mode  ipk,  where  the  energy 
is  in  the  sense  of  the  induced  norm  on  H.  Letting  P^pU  =  {u,(p)^  denote  the 
orthogonal  projection  of  u  onto  the  (normalized)  function  the  average  energy 
captured  by  9  is  given  by 

E{\\P^uf)  =  E{  {P^u,  P^u)  )  =  E{{  {u.  (u,  9>  ) 

=  E{{u,ip)  {ip,u)) 

=  E{\{u.^)\^), 

since  ||^||  =  1.  Thus,  by  (4,6),  Xk  is  the  average  energy'  captured  by  mode  ifk- 
Span  of  the  snapshots 

We  mentioned  in  the  previous  section  that  the  range  ol  R  —  E[u®u*)  is  contained 
within  the  span  of  the  ensemble  This  follows  directly  from  the  definition  of 

/?.  using  the  property  of  the  average  that  E{f{u))  is  a  linear  combination  of  /(u^). 
Note  that  in  the  stochastic  setting  which  is  often  used  for  discussing  the  POD,  the 
ensemble  of  u  is  infinite,  and  more  work  is  needed  to  show  this  [Holmes  et  ah, 
1996].  Note  also  that  the  range  is  not  necessarily  equal  to  the  span  of  the  as 
this  depends  on  the  average  (for  instance,  E{')  could  be  a  weighted  average,  giving 
zero  weight  to  certain  elements). 

An  important  imphcation  of  this  result  is  that  the  POD  modes  (clearly  in  the 
range  of  R)  may  be  written  as  linear  combinations  of  the  “snapshots”  That  is, 
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for  any  POD  mode  we  may  write 


m 

=  (4.10) 

fc=i 

for  some  scalar  coefficients  c^. 

Inherited  properties 

A  consequence  of  the  preceding  property  is  that,  if  every  snapshot  u'^  satisfies  a 
given  linear  constraint,  say  L{u'‘)  =  0  for  all  k,  where  L  :  H  R  is  linear,  then 
the  POD  modes  will  each  satisfy  that  constraint  as  well.  For  any  POD  mode  (p, 
by  (4.10)  we  have 


(m  \  m 

F  I  =  'F  =  0. 

k=i  /  k=i 

Examples  of  constraints  which  frequently  arise  are  the  divergence-free  condition 
for  velocit}^  in  incompressible  flow  (div  =  0)  and  finear  homogeneous  boundary 
conditions  (e.g,,  the  no  slip  condition  u*^(wall)  =  0). 

Symmetry 

In  the  presence  of  translational  s^mimetry  (i.e.,  if  the  data  is  homogeneous  in 
a  spatial  dimension),  it  can  be  shown  that  the  POD  modes  reduce  to  Fourier 
modes  in  the  directions  of  the  S3mimetry  [Sirovich,  1987].  Symmetries  do  not  arise 
in  the  flows  we  discuss  here,  so  we  do  not  dwell  on  this  point,  but  we  mention 
that  traveling  POD  modes  can  provide  more  accurate  results  (for  a  fixed  number 
of  mode*.':.)  than  the  standard  POD.  However,  if  reduced-order  models  are  to  be 
obtained  via  Galerkin  projection  (see  section  4.2),  then  one  needs  a  reconstruction 
equation  to  determine  the  position  of  the  traveling  wave,  and  to  close  the  resulting 
system  of  ODEs.  as  discussed  in  Rowley  L  Marsden  [2000]. 

4.1.3  Computation:  method  of  snapshots 

The  method  of  snapshots  provides  an  alternate  way  of  computing  the  POD  modes, 
and  is  often  more  efiicient  than  the  direct  method,  where  one  directly  solves  the 
eigenvalue  problem  (4.7). 

Let  the  average  E[^)  be  defined  as  a  weighted  average  over  the  snapshots 

£(/(»))  =  23 

k 

vhere  the  weights  ocj  >  0  satisfy  q ,  =  1.  (In  this  section,  all  sums  are  from 

i . T\T>ically,  Oj  =  1/m,  for  equal  weighting. 
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From  the  previous  section,  any  eigenfunction  ip  may  be  written  as  a  linear 
combination  of  snapshots 


ip^Y^Cku'^  (4.11) 

k 

for  some  coefficients  Ck  eC.  The  eigenvalue  problem  Rip  =  X(p  then  becomes 

=  A  ^  Cfcu'" 

k 

=  ^  AcfcU*' 
k 

=  ^Acfcu*. 
k 


/c  =  1, . . .  ,  m  (4-12) 

3 

which  is  just  the  standard  m-dimensional  eigenvalue  problem 

Uc=Xc,  (4.13) 

where  c  =  (cj, . . .  ,Cfn)  and  U  is  an  m  x  m  matrix  with  Uij  ~  ,u^).  Note 

that  U  is  self-adjoint  imder  the  inner  product  (a,  6)^  =  Ylj  eigen¬ 

functions  c  axe  orthogonal  under  the  same  inner  product. 

To  summarize,  the  direct  solution  of  (4.7)  involves  solving  an  eigenvalue  prob¬ 
lem  on  the  space  H,  which  may  be  infinite-dimensional.  For  instance,  if  the  snap¬ 
shots  e  H  are  data  from  a  finite-difference  simulation,  the  dimension  of  H  will 
be  the  number  of  gridpoints  in  the  simulation,  typically  large.  By  contrast,  the 
method  of  snapshots  involves  solving  an  m-dimensional  eigenvalue  problem  (4.13), 
where  m  is  the  number  of  snapshots  in  the  ensemble.  The  method  of  snapshots  is 
thus  more  efficient  whenever  the  number  of  snapshots  is  smaller  than  the  number 
of  gridpoints. 

4,1.4  Affine  subspaces 

The  POD  described  in  section  4.1  determines  a  /mear  subspace  5  which  optimally 
spans  an  ensemble  of  data.  Often,  it  makes  more  sense  to  project  the  data  onto 
an  affine  subspace.  Given  an  element  b  E  H  and  a  linear  subspace  5,  the  affine 
space  St  is  given  by 


E{u®u*)'^Cku'‘ 

k 

{u^  ®  {u^y)cku^ 

j  k 

It  suffices  to  solve 

^ak{u\u*‘)cj  =  Acfc, 


•56  =  {fc  +  I  S  5}.  (4.14) 

The  POD  procedure  extends  easily  to  affine  spaces,  and  amounts  to  subtracting 
the  element  b  from  the  ensemble  before  performing  the  POD. 
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Given  an  ensemble  of  data  {u^  €  i7},  and  an  element  6  G  i/,  we  form  the 
new  ensemble  where  u''  =  -  b.  The  standard  POD  then  provides  a  linear 

subspace  5  w'hich  optimally  spans  the  ensemble  {u*'}.  If  {ipj  |  j  =  1 . n}  is  the 

orthonormal  basis  for  S,  the  projection  onto  S  is,  as  before, 

n 

j=i 


and  the  projection  onto  5^  is  then 

n 

PstU  =  b  +  Y^{u-b,^j)ifj. 

3=\ 

Note  that  the  element  b  £  H  must  be  fixed  before  finding  the  POD  modes  . 
Tv'pically  we  will  choose  b  to  be  the  mean  of  the  ensemble,  so  b  =  E{u).  but  other 
choices  are  possible. 

4.2  Galerkin  projection 

Consider  a  djmamical  system  which  evolves  in  a  Hilbert  space  H.  In  particular, 
for  u{t)  €  H.  u{t)  satisfies 


■u  =  X{u).  (4.15) 

where  A  is  a  vector  field  on  H.  For  instance,  for  a  partial  differential  equation 
(PDE)  governing  a  variable  u{x.t),  defined  on  some  spatial  domain  x  €  fi,  will 
be  a  space  of  functions  defined  on  D,  and  A'  will  be  a  spatial  differential  operator. 
Given  a  finite-dimensional  subspace  S  of  H,  we  wish  to  determine  a  dynamical 
system  which  evolves  on  S  and  approximates  (4.15)  in  some  sense.  This  dynamical 
system  will  be  denoted 


r  =  Xs{r),  (4.16) 

where  r{t)  C  S  and  A5  is  a  vector  field  on  S.  Galerkin  projection  specifies  this 
vector  field  by 


A'5(r)  =  PA'(r).  (4.17) 

w^here  F  :  H  —*  S  is  the  orthogonal  projection  map.  The  projection  theorem 
states  that  this  approximation  is  optimal,  in  the  sense  that  it  minimizes  the  error 
||A5(r)  —  A(r)||,  where  j|  •  ||  is  the  induced  norm  on  H  [Luenberger,  1997]. 

Galerkin  projection  is  a  special  case  of  weighted  residual  methods  [Fletcher, 
1991],  which  specify  the  vector  field  Xs  by  requiring 


{Xs{t)  -  A'(r),u'fc>  =  0. 


fc  =  1, . . .  ,n. 


(4.18) 
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for  some  specified  weight  functions  G  i/,  where  n  =  dim  5..  This  is  a  linear 
system  of  n  linear  equations  which  determine  the  n-dimensional  vector  Xs(r).  In 
the  Galerkin  method,  the  weight  functions  Wk  are  chosen  to  be  basis  functions  for 
the  subspace  5,  and  so  the  residual  Xs(r)  -  X(r)  is  orthogonal  to  5,  which  is 
equivalent  to  (4.17). 

To  apply  this  method  to  computations,  we  need  to  write  (4.16)  in  coordinates. 
Let  {ipk  €  1  A:  =  1, ...  ,72}  be  a  basis  for  a  subspace  5.  (For  instance,  such  a 

basis  may  be  obtained  from  POD  on  a  set  of  data,  as  described  in  the  preceding 
section.)  Writing  r{t)  in  coordinates  aj^{t)  with  respect  to  this  basis,  w^e  have 

n 

r{t)  =  '^akit)ipk,  (4.19) 

k=l 

and  using  (4.18)  with  (4.16)  becomes 

n 

'^dj{t){<pj,<pk)  =  {X{r{t)),<pk) ,  (4.20) 

j=l 

This  defines  a  hnear  system  of  ODEs  for  the  evolution  of  aj{t).  If  the  basis  func¬ 
tions  are  orthonormal,  then  (4.20)  further  simplifies  to 

dk{t)  =  {X{r{t)),pk) .  A—  1,...  ,72.  (4.21) 

This  is  why  an  orthonormal  basis  is  especially  desirable,  because  we  avoid  inverting 
an  72  X  n  matrix  to  solve  for  a. 

The  Galerkin  procedure  also  extends  easily  to  affine  spaces,  mentioned  in  sec¬ 
tion  4.1.4.  Given  the  linear  subspace  5  as  above,  and  an  element  6  G  i/,  we  wish 
to  project  the  dynamics  onto  the  affine  space  ^  {b  v  \  v  £  S}.  In  this  case, 
the  expansion  (4.19)  takes  the  form 

n 

r{t)  =  6  +  ^  ak{t)ipk-  (4.22) 

A;=l 

where  r{t)  G  5^.  and  the  projected  equations  are  the  same  as  above,  (4.20) 
or  (4.21).  the  only  difference  being  the  new  expression  (4.22)  for  r{t). 

4.2.1  Example:  Quadratic  equations 

For  many  t}q)es  of  equations,  the  ODEs  given  by  (4.21)  may  be  determined  analyt¬ 
ically,  in  terms  of  the  coordinates  This  is  particularly  useful  computationally, 
as  the  inner  product  in  (4.21)  does  not  need  to  be  computed  at  every  timestep. 
We  illustrate  this  with  a  quadratic  PDE  described  by  u  =  X{u),  with 

X{u)  —  L{u)  -h  Q(u,  u). 
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where  L  :  H  H  is  a.  Hneax  operator,  and  Q  :  H  x  H  H  is  a  bilinear  operator 
(linear  in  each  argument).  W'e  project  the  dynamics  onto  the  aiRne  space  Sb, 
expanding  r{t)  €  as  in  (4.22).  The  projected  ODEs  (4.21)  become 

dk{t)  =  {L{r)  +Q{r,r),(pk) 

L(b  +  Y^a,{t)Lpi'^  +  + 

•  i  j 

—  +  y^Q]fa,(t)aj(t), 

where  the  quantities 

h  =  {L{b)  +  Q{b,b),^k) 

Lk  =  m^i)  +  Q{b,ipi)  +  Q(<fi,b),ipk) 

Qk  =  iQi‘Pi’‘Pj)^Vk) 

are  constants  (independent  of  t)  which  may  be  determined  before  integrating  the 
ODEs. 

A  similar  simphfication  is,  of  course,  possible  for  cubic  and  higher-order  poly¬ 
nomial  nonlinearities,  but  for  higher  powers,  the  number  of  coefficients  becomes 
unwieldy  for  even  moderate  dimensions  n,  so  we  restrict  ourselves  to  quadratic 
governing  equations. 

4.3  Incompressible  fluids 

The  applications  in  this  thesis  concern  compressible  Sows,  but  we  outline  the 
well-established  methods  for  applying  POD/Galerkin  to  incompressible  flows,  and 
highlight  the  differences  for  compressible  flows  in  the  next  section. 

The  definition  of  incompressibility  is  that  the  velocity  u  =  {u.v,w)  is  diver¬ 
gence  free,  div  u  =  0.  The  motion  of  the  fluid  satisfies  the  Navier-Stokes  equations 
^Batchelor.  1967].  written  in  Cartesian  coordinates  as 

Du  _  _o 

—  =  -Vp  -!-  uV'u,  (4.23) 

where  if  is  the  viscosity,  p  is  the  pressure,  and  D/Dt  =  d/dt  +  u  •  V  is  the  material 
deri\'ative.  Velocities  have  been  normalized  by  some  velocity  scale  17,  lengths  by  a 
length  scale  I,  time  by  U/L,  pressure  by  pU^  where  p  is  the  density,  and  viscosity 
by  pUL  (thus,  u  is  the  reciprocal  of  the  Reynolds  number).  These  equations  may 
be  written  as 

ii  =  A''(u)  -  Vp 


where  A'(u)  =  -(u  •  V)u  -|-  t/V^u. 


(4.24) 
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For  incompressible  flows,  we  write  the  velocity  u  as  an  expansion  in  POD 
modes  9?(x),  defined  on  a  spatial  domain  fl  in  which  the  fluid  evolves: 

n 

=  (4-25) 

j=l 

Our  Kilbert  space  H  is  just  the  space  of  smooth  (C^)  vector- valued  functions 
on  n,  with  the  standeird  inner  product 

(u,v)=  [  u(x)*v(x)dK  (4.26) 

Jq 

(Additionally,  we  restrict  the  space  H  to  contain  only  functions  with  finite  norm.) 

Inserting  the  expansion  (4.25)  into  the  Navier-Stokes  equations  (4.24),  and 
taking  an  inner  product  with  gives 


dk  =  {N{u),  -  (Vp,  (pk) .  (4.27) 

The  pressure  term  on  the  right-hand  side  is  easily  rewritten 

(^P-^k)=  /  <Pk-'^pdV=  /  div{pipk)dV  =  I  pipi^-ndS 

Jq  Jq  JdQ 

where  in  the  second  equality  we  have  used  that  divc^^.  =  0.  Thus,  this  term 
depends  only  on  the  pressure  on  the  boundary  dQ.  Furthermore,  if  velocity  is  zero 
along  the  boundary  (for  instance,  at  a  wall,  or  in  the  farfield  of  an  open  flow), 
then  =  0  on  dfl,  and  the  pressure  term  vanishes  altogether.  If  the  boundary  is 
not  a  wall,  but  rather  an  ‘‘artificial”  boundetry  we  impose  (i.e.,  we  consider  only  a 
limited  po  iion  of  the  whole  flow'),  then  the  pressure  term  represents  the  influence 
of  the  rest"  of. the  flow’  on;  the  domain  we  are  considering,  and  must  be  specified  as 
a  boundary  condition. 

The  technique  for  incompressible  flow^s  is  w^ell  know’n.  Many  enhancements  to 
the  basic  theory  exist  (see  Holmes  et  al.  [1996]  for  a  thorough  discussion),  but  the 
essential  elements  are  the  same  as  given  here.  The  main  feature  that  is  normally 
included,  but  is  not  mentioned  above,  is  that  the  velocity  is  usually  decomposed 
into  mean  and  fluctuating  components  (u  =  u  +  u'),  where  the  mean  u  is  slowdy 
varying  in  time.  The  mean  may  then  be  described  in  terms  of  the  fluctuations  u', 
w’hich  give  rise  to  Reynolds  stresses,  as  discussed  in  Aubry  et  al.  [1988].  This  is  an 
important  distinction,  and  gives  rise  to  cubic  equations,  w’hile  the  equations  (4.27) 
are  only  quadratic.  Another  enhancement  to  the  basic  theory  is  the  modeling  of 
the  energy'  transfer  to  the  higher  modes.  In  this  thesis,  we  take  the  mean  flow  to  be 
constant  in  time,  w’ithout  attempting  to  model  the  Reynolds  stress  contributions, 
and  w’e  neglect  the  energ\'  transfer  to  the  higher  modes.  These  extensions  could 
presumably  be  added  to  the  compressible  theory  as  well. 
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4.4  Compressible  fluids 


One  of  the  main  contributions  of  this  thesis  is  a  method  for  applying  the  above 
techniques  to  compressible  flows.  The  distinction  is  a  fundamental  one.  On  a 
superficial  level,  the  constraint  divu  =  0  no  longer  holds.  On  a  more  fundamental 
level,  in  a  compressible  flow,  the  thermodynamic  variables  are  dyamically  impor¬ 
tant,  and  must  be  included  in  the  actual  dynamics,  not  merely  as  a  constraint. 
For  incompressible  flows,  as  discussed  above,  the  pressure  drops  out  completely,  or 
appears  only  as  an  imposed  boundary  condition,  and  the  only  dynamical  variable 
is  the  velocity.  For  compressible  flows,  this  is  not  the  case,  and  evolution  equa¬ 
tions  must  be  given  for  one  or  more  thermodynamic  variables  (e.g.,  p,  p,  entropy  s, 
enthalpy  /i),  as  well  as  the  velocity.  This  introduces  further  questions  of  whether 
to  treat  these  variables  completely  separately  from  the  velocity,  or  together  as  a 
single  vector-valued  “configuration”  variable  (e.g.,  q(x)  =  (p,u,z;,p)(x)). 

We  consider  both  alternatives,  first  discussing  a  scalar-valued  method,  where 
we  compute  separate  POD  modes  for  each  flow'  variable,  and  next  discussing  a 
vector-valued  method,  where  we  write  all  the  flow  variables  together  as  a  .single 
vector,  and  compute  one  set  of  vector-valued  POD  modes.  To  do  this,  we  first  need 
to  define  an  appropriate  inner  product  on  our  configuration  space — the  standard 
inner  product  may  not  be  a  sensible  choice.  For  instance,  if  we  use  as  flow  variables 
q  =  (p,  u.  t',p),  defined  on  the  fluid  domain  fl.  the  standard  inner  product  is 


(qiiq2)=  /(PiP2  +  UiU2  + 


(4.28) 


which  does  not  make  dimensional  sense,  since  one  cannot  add  a  velocity  ^md  a 
pressure.  We  introduce  an  inner  product  suitable  for  compressible  flows  in  sec¬ 
tion  4.4.3. 


4.4.1  Governing  equations 

The  compressible  Navier  Stokes  equations,  written  in  Cartesian  tensor  notation, 
are  given  by 


Dp 

Dt 

Du, 

'W 

Dt 


-f  p  div  u  =  0 


+  ^2p  5, 


dp 

dxj  dxj‘ 


-Sij  divu 


1  /  dui  duj  \ 

2\^j  eir,) 


(see  Batchelor  [1967]),  where 
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is  the  rate-of-strain  tensor  (the  symmetric  part  of  the  velocity  gradient),  and 

is  the  rate  of  dissipation  of  mechanical  energy.  Here,  s  is  the  entrop}^  T  is  the 
temperature,  and  k  the  thermal  conductivity  of  the  fluid. 

These  equations  are  significantly  more  complicated  than  the  incompressible 
equations,  so  we  make  some  approximations  to  obtain  a  simpler  set  of  equations 
which  we  may  use  for  Galerkin  projection.  For  the  cavity,  we  consider  cold  flow 
(^wall  =  Too)  and  note  that  if  the  Mach  number  is  not  too  high,  density  gradients 
will  remain  small  and  will  be  dominated  by  presure  changes.  This  is  consistent  with 
the  neglect  of  the  viscous  dissipation  $  and  heat  conduction  in  the  energ>^  equation, 
and  thus  we  treat  the  flow  as  isentropic.  However,  we  retain  the  viscous  diflFusion 
in  the  momentum  equation,  but  in  this  term  we  assume  dynamic  and  kinematic 
viscosities  are  constant,  again  under  the  approximation  that  temperature  gradients 
are  small.  Under  these  assumptions,  the  equations  of  motion  become 


where  i/  =  is  a  constant.  We  still  have  more  variables  than  equations,  so  to 
close  the  system  we  require  an  equation  of  state.  With  an  equation  of  state,  we 
may  write  the  equations  in  terms  of  a  single  thermod\mamic  variable,  since  the 
entropy  is  constant.  Any  thermod^mamic  variable  (independent  of  entropy)  will 
suffice,  but  the  enthalpy  h  is  particularly  convenient,  as  we  shall  see  shortly. 
Using  the  ideal  gas  relation  p  =  pRT,  w’e  obtain  the  relation 

=  (^  _  l)h,  (4.30) 

P 

where  q  is  the  local  sound  speed.  Using  the  Gibbs  equation 

dh  =  Tds  +  -dp  (4.31) 

P 

with  ds  =  0,  the  equations  (4.29)  are  easily  written  in  terms  of  h  alone: 


-h  (")  -  l)/2divu  =  0 


Du 

Dt 


-f-  V/i  =  i^V^u. 


(4.32) 


This  procedure  will  work  for  any  thermodynamic  variable,  but  the  enthalpy  is 
particularly  convenient,  because  the  equations  (4.32)  are  quadratic.  If  we  used  the 
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density,  for  instance,  the  momentum  equation  would  become 


-^  +  p'’’  ^Vp=i/V^u. 

For  air,  7  =  1.4,  so  the  exponent  7  —  2  in  the  above  equation  is  fractional,  which 
would  cause  difficulty  when  we  perform  Galerkin  projections,  described  in  sec¬ 
tion  4.2. 

Another  choice  of  thermodynamic  variable  which  yields  quadratic  equations  is 
the  local  sound  speed  a.  Using  (4.30),  the  equations  (4.32)  become 


Da  7  -  1 
Du  2 


a  div  u  =  0 
aVa  =  i/V^u. 


(4.33) 


The  equations  (4.32)  and  (4.33)  are  equivalent,  and  we  refer  to  them  (in  either 
form)  as  the  isentropic  Navier-Stokes  equations.  We  use  the  enthalpy  form  (4.32) 
with  the  scalar-valued  method  discussed  in  the  next  section,  and  the  equations 
in  terms  of  sound  speed  (4.33)  with  the  vector-valued  method  discussed  in  sec¬ 
tion  4.4.4. 


4.4.2  Scalar-valued  POD  modes 

In  this  method,  we  consider  each  flow  variable  separately,  and  expand  each  flow 
variable  in  its  own  set  of  POD  modes: 


u{x,y.t)  =  u{x.y)  +  'Y^Uj(t)id  (x.y) 
j=l 

v{x.y,t)  =  v(x.y)  +  J2vj(t)vJ(x.y)  (4.34) 

J=l 

n 

h{x,y.t)  =  h{x.y)  +  J2kj{t)M{x.y), 
j=i 

where,  c.g..  u^ix.  y)  are  the  POD  modes.  Uj{t)  are  the  time  coefficients,  and  u(x,  y) 
is  fixed,  usually  taken  to  be  the  mean  of  the  snapshots  used  for  the  POD.  The  inner 
product  we  use  here  is  just  the  standard  inner  product  on  real-valued  functions: 

if -9)-  [  f{x.y)g(x,y)dV 
JQ 

where  Q  C  is  the  region  in  which  the  fluid  evolves. 

.Next,  wo  write  the  isentropic  Navier-Stokes  equations  (4.32)  in  nondimensional 
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form,  as 

Ut  -f  UUx  +  VUy  +  hx  =  (“Uxx  +  ^yy) 

Vt  -f  UVj:  +  VVy  hy  =  +  Vyy) 

ht  +  uha:  +  vhy  +  (7  ~  l)h{Ux  +  Vy)  0, 


(4.35) 


where  Re^  =  VLjv.  Here,  velocities  have  been  normalized  by  the  freestream 
velocity  U,  h  has  been  nondimensionalized  by  a:  and  y  by  a  length  L,  and  time 
by  L/U.  Inserting  the  expansions  (4.34)  into  these  equations,  and  taking  inner 
products  with  the  POD  modes  gives  the  system  of  ODEs 

n  n 

^k  =  bk  +  Yl  +  ‘^kVi  +  iVijkUiUj  +  q'tjkViUj) 

1=1  ij=l 

n  n 

^’k  =  bl  +  ^  (4ui  +  (Fik'^i  +  e\^hi)  +  ^  {v'tjk'^iVj  +  q'ijk'^iVj)  (4.36) 

i=l  i,j=l 

n  n 

bk  —  b*^  +  {Pijk'^ibj  +  q'^j^.Vihj), 

t=i  i,i=i 

where  the  coefRcients  are  given  by 

bk  =  -  (uux  +  VUy  +  hx-  {uxx  +  Uyy)/ReL, 
b^  —  ^UVx  “1“  by  {_Vxx  4"  ^’yy)/R^jLi  ^  ) 

fcjt  =  -  (uhx  +  vhy  +  (7  -  l)h{ux  +  Vy),  h':  'j  . 

cDt  =  -  {uu\  +  vil'y  +  u'Ux  -  +  Wyy)/^^/,,  U*") 

cu- =  - (w’Ci, €■*■■) 

c':,  =  -{u%+{^-i)hu)x.h>^) 

<4  =  - 

<^*r  =  -  +  vvl  +  v'vy  -  +  vly)/ReL,v^'^ 

d^k--{v%+h-l)H,b'^) 


70 


4  Model  Reduction  for  Fluids 


er, u'^-) 

^ik  =  -  +  vhl  +  (7  -  +  Vy), 

Pxjk  =  {u'ui,u'^)  q^jk  = 

Ptjk  =  v’')  q'ljk  =  v’') 

pU  =  +  (7  -  l)h^ui,h^)  q^k  =  +  (7  -  l)h^vi,  . 

The  equations  (4.36)  represent  a  reduced-order  model  of  the  isentropic  Navier- 
Stokes  equations.  The  original  system  (4.32),  consisting  of  three  coupled  PDEs, 
has  been  reduced  to  a  system  of  3n  coupled,  quadratic  ODEs,  where  n  is  the 
number  of  POD  modes  used  for  each  variable.  For  many  situations,  the  number  of 
modes  n  required  to  capture  a  large  fraction  of  the  energy  is  small  (say  <  20),  so 
these  equations  represent  a  significant  reduction  in  order,  since  the  typical  number 
of  gridpoints  used  in  a  simulation  is  0(10®).  We  apply  these  equations  to  the  flow 
over  a  rectangular  cavity  in  chapter  5. 

Note  that  we  could  leave  Rej;,  as  a  parameter  in  the  equations  (4.36),  but  in 
this  thesis  we  do  not  study  the  bifurcation  behavior  of  the  equations  for  scalar¬ 
valued  POD  modes.  We  do  include  the  effect  of  parameters  for  the  equations  for 
vector-valued  modes,  however,  discussed  in  section  4.4.4. 

4.4.3  Inner  products  for  compressible  flow 

In  order  to  obtain  vector-valued  POD  modes,  and  perform  Galerkin  projections  as 
in  the  previous  section,  we  must  first  define  an  inner  product.  Here  we  introduce  a 
family  of  inner  products  useful  for  compressible  flow  problems.  As  mentioned  ear¬ 
lier.  the  usual  inner  product  (e.g.,  (4.28))  may  not  make  dimensional  sense  when 
both  thermodynamic  and  kinematic  variables  are  included.  Of  course,  one  could 
simply  nondimensionalize  the  variables,  but  then  the  sense  in  which  projections 
are  "optimal  is  rather  arbitrary,  and  depends  on  the  nondimensionalization  (es¬ 
sentially  the  choice  of  "coordinates'’),  not  on  any  physical  grounds.  We  seek  an 
inner  product  which  makes  intuiti%-c  sense,  in  that  the  "energy”  defined  by  the 
induced  norm  is  a  meaningful  physical  quantity. 

It  is  natural  for  the  induced  norm  to  be  related  to  the  total  energy  of  a  system. 
For  instance,  for  a  mechanical  system  this  might  be  the  kinetic  energy — indeed,  for 
incompressible  flows,  the  norm  induced  by  the  stanard  inner  product  is  the  kinetic 
energy.  For  compressible  flows,  both  thermodynamic  and  kinematic  variables  con¬ 
tribute  to  the  total  energy.  For  instance,  the  stagnation  (or  total)  enthalpy  of  the 
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flow  is  given  by 

ho  = 

where  h  is  the  static  enthalpy.  Analagously,  the  stagnation  energy  is  e  =  £"  -r 
where  E  =  h/'y  is  the  internal  energy  per  unit  mass.  Motivated  by 
these  physical  quantities  based  on  energy,  we  look  for  inner  products  which  ha\^ 
induced  norms  of  the  form 

^llqfa  =  +  (4.37) 


where  q  is  the  vector  of  flow  variables,  and  a  >  0  is  a  constant.  Note  that  the  right- 
hand  side  of  (4.37)  is  not  quadratic,  because  h  appears  linearly,  but  this  is  easily 
remedied  by  transforming  to  the  flow  variables  q  =  since  —  (7  —  l)/i, 

by  (4.30).  Thus,  we  define  a  family  of  inner  products 


(qi>q2)< 


-IJ 


U1U2  -f  V1V2  + 


2a 

7  -  1 


aia2 )  dV, 


(4.38) 


which  has  induced  norm  given  by  (4.37). 

Obviously,  choosing  a  =  1  corresponds  to  using  the  integral  of  the  stagnation 
enthalpy  as  the  norm,  and  this  is  the  choice  we  take  for  the  results  presented  in 
this  thesis.  Taking  a  =  I/7  corresponds  to  using  the  stagnation  energy  as  the 
norm.  However,  note  that  neither  the  integral  of  the  stagnation  enthalpy  or  the 
stagnation  energ>’  is  actually  a  conserved  quantity.  The  conserved  quantity  is  the 
total  energy,  given  by 


pE  +  ^p{u^  +  dV:  ..(4.39) 

(Of  course,  this  is  only  conserved  if  there  is  no  energy  flux  through  the  boundary 
of  fl,  as  when  the  velocity  goes  to  zero  on  9Q.)  It  is  not  obvious  how  to  define  an 
inner  product  which  has  (4.39)  as  the  induced  norm,  though  this  would  perhaps 
be  the  most  natural,  from  an  energ>'  point  of  view.  However,  we  note  that  for 
Q  =  l/(")  -  1),  the  induced  norm  actually  is  a  conserved  quantity,  for  inviscid, 
isentropic  flows. 

With  Q  —  1/(7  —  1),  the  integrand  of  (4.37)  becomes 

C=— ft  +  ju.u, 

where  u  =  (u.u).  Thus,  assuming  isentropic  flow  (equations  (4.32)  with  t/  =  0), 
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we  have 

DC  _  1  Dh  Du 

Dt 

=  “/idivu  —  u  ■  V/i 
=  — div(/iu). 

Integrating  over  Q  and  using  the  divergence  theorem  then  gives 

Hence,  if  u  =  0  along  50,  then  the  integral  is  zero,  and  the  induced  norm  is 
conserved. 

4.4.4  Vector- valued  POD  modes 

We  begin  by  nondimensionalizing  the  equations  (4.30),  scahng  the  velocities  u,v 
b\  the  freestream  velocity  C/,  the  local  sound  speed  by  the  ambient  sound  speed 
X  and  y  by  a  length  L,  and  time  by  L/U.  The  equations  become 

Ut  +  UUx  +  VUy  +  +  Uyy) 

Vi  +  uv^  +  vvy  +  J^^auy  =  ^^(Uxi  +  Vyy)  (4.40) 

at  +  uox  +  VUy  +  ^a{ux  +  Vy)  =  0, 

where  Re^  =  UL/v  and  M  —  U/a^o-  This  nondimensionalization  is  particularly 
useful,  because  the  Mach  number  appears  in  the  equations.  Thus,  M  will  also 
appear  in  Galerkin  projectioas,  and  we  may  investigate  its  effect  on  the  resulting 
reduced-order  models.  Note  that  by  scaling  x  and  y  by  different  quantities  (e.g., 
scale  X  by  cavity  length  L  and  y  by  momentum  thickness  0),  we  introduce  another 
nondimensional  parameter,  the  ratio  of  these  lengths.  For  the  cavity  flow,  this 
procedure  may  be  used  to  investigate  effects  of  L/Oq,  for  instance,  but  in  this 

thesis  we  consider  only  the  effects  of  M  and  Re/..  Writing  q  =  (u,t;,a),  the 

equations  may  be  written 

^  =  ^^(q)  +  ](^<?i(q-q)  +  <32(q.q)  (4.41) 

where 

=  (4,42) 


VXX  "b  Vyy  1  , 
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Q2(q\q^)  =  - 


/v}u\  + 

I  U^v'l  +  V'^Vy 
yu^a^j.  +  v^Oy  + 


^a^ul 


^Ih 


Next,  we  expand  q  in  terms  of  POD  modes,  as 


(4.43) 


q(x,t)  =  q(x)  +  ^aj(f)v5j(x),  (4.44) 

j=i 


where  q  is  fixed,  and  typically  taken  to  be  the  mean  of  all  the  snapshots  used 
for  POD.  Recall  that  q  must  be  subtracted  from  the  original  data  before  the 
POD  modes  are  computed,  as  described  in  section  4.1.4.  The  resulting  Galerkin 
equations  are  given  (see  section  4.2)  by 


1 

i,j=l 


where  the  coefficients 


(4.45) 


6'  =  (L(q),<^,) 

=  (Qi(q-q)-¥’fc> 
tfc  =  (Q2(q-q)-¥’i) 


Lj,  =  (L(<^.).<^,) 

=  (<?i(q-v5.)  +  Qii‘Pf^)^‘Pk) 

L%  =  (<32(q.  <Pi)  +  Q2{(pi,q)- ¥>k) , 


Qljk  =  {Qiiv>i^v>j)^‘Pk) 

Q^jk  =  {Q2i‘Po‘Pj)<<pk) 

are  constants  which  may  be  computed  before  solving  the  reduced  system.  We  use 
the  inner  product  from  section  (4.4.3)  (with  a  —  1)  both  for  computing  the  POD 
modes  and  for  the  Galerkin  projection. 

The  equations  (4.45)  represent  the  reduced-order  model  of  the  isentropic  Navier- 
Stokes  equations  when  vector-valued  modes  are  used.  They  are  analogous  to  equa¬ 
tions  (4.36)  for  scalar  modes.  Note  that  (4.45)  is  a  system  of  n  equations,  while 
(4.36)  is  a  system  of  3n  equations,  where  n  is  the  number  of  modes  used  in  the 
Galerkin  projection.  Thus,  it  is  reasonable  to  expect  that  lower-order  models  will 
be  possible  with  (4.45)  than  with  (4.36).  For  the  cavity  flow,  this  is  indeed  the 
case,  as  we  shall  see  in  chapter  5. 
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Chapter  5 

Low-order  Models  of  Cavity  Oscillations 


The  model-reduction  techniques  discussed  in  the  previous  chapter  are  particularly 
well  suited  to  cavity  flows.  First,  the  dynamical  behavior  of  cavity  flows  is  relatively 
simple,  and  exhibits  qualitative  features  of  low-dimensional  dynamical  systems 
(e.g.,  bifurcations,  limit  cycles),  so  it  is  reasonable  to  expect  that  a  low-dimensional 
description  is  possible.  Second,  as  mentioned  in  the  introduction,  one  of  the  main 
reasons  one  is  interested  in  modeling  is  to  better  understand  how  to  control  a 
flow.  In  many  situations,  the  goal  of  control  is  to  drive  the  flow  into  a  completely 
different  state — for  instance,  relaminarization  of  a  turbulent  flow — ^which  no  longer 
resembles  the  original,  uncontrolled  flow.  In  these  situations,  the  POD/Galerkin 
method  may  not  be  suitable,  since  the  POD  modes  capture  only  the  limited  portion 
of  phase  space  that  is  spanned  by  the  snapshots.  For  the  cavity  flow,  however,  we 
expect  the  controlled  flow  to  have  the  same  general  flow  features  of  the  uncontrolled 
flow:  a  shear  layer  will  still  span  the  cavity,  and  will  still  amplify  disturbances, 
but  the  amplitude  of  the  oscillations  will  be  smaller,  or  ideally,  will  vanish. 

Th^'  ^oal  of  this  chapter  is  to  apply  the  POD/Galerkin  method  to  the  data  from 
our  simulations,  to  develop  reduced-order  models  for  the  cavity  flow^  These  models 
are  suitable  for  bifurcation  analysis,  but  do  not  yet  include  effects  of  actuation,  and 
so  are  not  yet  directly  suitable  for  control  design  or  analysis.  We  use  both  scalar¬ 
valued  and  vector-valued  POD  modes,  and  show  that  the  vector-valued  method 
has  significant  advantages  over  the  scalar- valued  method.  Throughout,  w'e  focus 
on  two  runs  from  table  3.1:  runs  L2  and  H2.  Both  of  these  runs  are  oscillating 
in  shear-layer  mode,  but  run  L2  has  oscillations  at  both  of  the  first  two  Rossiter 
frequencies,  while  run  H2  has  oscillations  only  at  Rossiter  Mode  II. 


5.1  Scalar- valued  modes 

All  of  the  results  in  this  section  pertain  to  the  methods  for  scalar-valued  POD 
modes,  discussed  in  section  4.4.2. 
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Run 

L/D 

L/9 

Ree 

istart 

lend 

iskip 

L2 

2 

52.8 

56.8 

20,750 

44,000 

250 

H2 

2 

58.4 

68.5 

70,000 

80,000 

200 

TK4b 

4 

60.2 

58.8 

30,000 

45,000 

200 

Table  5.1:  Parameters  for  different  runs  considered,  along  with  the  timesteps  cf  the 
first  and  last  snapshots  (istart,  lend)  and  number  of  timesteps  between  snapshots 
(iskip). 


5.1.1  POD  modes 

We  compute  POD  modes  for  each  flow  variable  separately,  according  to  the  ex¬ 
pansion 


q{x,  y,  t)  =  q{x,  y)  +  ^  qk{t)qk{x,  y)  (5.1) 

A:=l 

where  q  is  any  one  of  {u, Qkix.y)  are  the  POD  modes,  qk{t)  the  time 
coefficients,  and  the  mean  q  is  taken  to  be  the  time  average  of  the  snapshots.  The 
method  of  snapshots  is  used,  and  the  matrix  U  from  equation  (4.13)  is  saved  so 
that  it  need  not  be  recomputed  if  one  wishes  to  enlarge  the  ensemble  of  snapshots. 

For  the  computation  of  scalar- valued  POD  modes  in  this  section,  snapshots  are 
first  interpolated  onto  a  uniform  grid  {xi,yj),  with  grid  spacing  (Ax,  A?/),  and  the 
inner  product  used  is  the  standard  Euclidean  inner  product,  given  by 

if^g)  =  '^f{xi.yj)g{xr,yj)AxAy,  (5.2^ 

for  any  functions  f^g.  The  uniform  grid  is  coeirser  than  the  original  grid  in  the 
DXS,  and  the  domain  is  truncated.  For  runs  wdth  L/D  —  2,  the  coarse  grid  has 
101  X  51  gridpoints  inside  the  cavity  (for  (x,y)  e  [0,2]  x  [-1,0]),  and  201  x  101 
gridpoints  above  the  cavity  (for  (x.y)  €  [-1.3]  x  [0,2]).  Spatial  derivatives  are 
computed  on  the  fine  grid  and  then  interpolated  onto  the  coarse  grid,  for  better 
accuracy  in  the  Galerkin  equations. 

The  snapshots  used  for  each  set  of  POD  modes  are  listed  in  table  5.1,  along 
with  the  relevant  parameter  values.  We  also  compute  POD  modes  for  a  third  run, 
listed  in  table  3.1  as  TK4b,  which  will  be  used  later  to  investigate  the  relative 
importance  of  the  parameters  L/D  and  L/6q.  Figure  5.1  shows  the  u- velocity  at  a 
point  in  the  shear  layer,  and  indicates  the  times  w^here  the  snapshots  were  taken. 
Snapshots  for  run  L2  included  more  of  the  developing  region,  as  shown. 

Run  L2 

Figure  5.2  shows  the  eigenvalues  for  the  POD  modes  for  run  L2.  Recall  that  the 
eigenvalue  X^  represents  the  ^‘energ>-’  captured  by  POD  mode  k,  in  the  sense  of  the 
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Figure  5.1:  Normal  velocity  at  a  point  in  the  shear  layer,  for  each  of  the  three 
runs:  (o)  indicates  where  snapshots  are  taken. 


Figure  5.2:  Eigenvalues  from  POD  of  run  L2:  u  (o),  v  (□),  h  (0),  and  lj  (a) 
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Mode  # 

%  energy  {u) 

%  energy  (i;) 

%  energy  {h) 

%  energy  lj 

1 

20.85 

21.39 

25.64 

20.37 

2 

19.34 

17.92 

21.28 

13.33 

3 

9.08 

10.29 

10.06 

8.19 

4 

8.41 

9.59 

9.24 

7.22 

Table  5.2;  Fraction  of  energy  captured  by  each  POD  mode,  run  L2. 

induced  norm  ||  •  |p  =  (•,  •),  which  here  is  given  by  (5.2),  where  each  flow  variable 
(u,  V,  h.  u)  has  its  own  norm.  Note  that  the  vorticity  has  significantly  more  energ}'. 
since  differentiation  amphfies  the  small  scales.  The  velocities  u.  v  have  more  energy 
than  the  thermodynamic  variable  h,  and  all  variables  fall  off  fairly  rapidly,  with 
the  energ>’  decreasing  by  an  order  of  magnitude  after  only  10  modes.  The  fraction 
of  the  total  energy  captured  by  each  of  the  first  4  modes  is  indicated  in  table  5.2. 

The  mean  of  the  snapshots  is  shown  in  figure  5.3,  and.  shows  the  shear  layer 
spanning  the  cavity,  with  a  steady  vortex  in  the  downstream  half  of  the  cavity. 

The  first  4  POD  modes  are  shown  in  figures  5.4-5. 7.  The  separate  POD  modes 
for  the  individual  variables  all  have  similar  features:  Modes  1  and  2  have  slightly 
less  than  one  full  wavelength  across  the  shear  layer,  corresponding  to  Rossiter 
Mode  I  (see  section  3.2.4),  and  modes  3  and  4  have  slightly  less  than  two  wave¬ 
lengths  across  the  shear  layer,  corresponding  to  Rossiter  Mode  II.  This  correspon¬ 
dence  is  also  confirmed  by  the  time  coefficients  q{t)  in  (5.1):  the  time  coefficients 
are  approximately  sinusoidal,  with  the  first  two  coefficients  at  the  first  Rossiter 
frequency,  and  the  next  two  coefficients  at  the  second  Rossiter  frequency.  Note 
also  that  the  POD  modes  arise  in  pairs,  out  of  phase  by  90°,  much  like  sine  and 
cosine,  with  Modes  1  and  2  forming  one  pair,  and  Modes  3  and  4  forming  another. 
(.Note  that  the  sign  of  the  POD  modes  is  arbitrary.) 

Run  H2 

The  POD  modes  for  run  H2  are  computed  the  same  way,  and  the  eigenvalues  are 
shown  in  figure  5.8.  The  energ>'  falls  off  much  more  steeply,  dropping  an  order  of 
magnitude  after  the  first  four  POD  modes.  This  faster  dropoff  is  to  be  expected, 
since  run  L2  has  a  more  complicated  behavior.  Also  note  the  pairwise  dropoff  of 
the  energy.  The  percent  energy  captured  in  each  mode  for  the  first  four  modes  is 
given  in  table  5.3. 

The  mean  and  first  four  POD  modes  are  plotted  in  figures  5.9-5.13.  From  the 
mean,  it  is  apparent  that  the  steady  vortex  in  the  rear  of  the  cavity  is  significantlv 
stronger  than  in  run  L2.  Also,  by  contrast  with  run  L2,  POD  modes  1  and  2  have 
nearly  two  full  wavelengths  in  the  shear  layer,  corresponding  to  Rossiter  Mode  II. 
Again,  this  correspondence  is  confirmed  by  the  frequency  of  the  time  coefficients. 
POD  modes  3  and  4  contain  smaller  spatial  scales,  possibly  corresponding  to  a 
higher  Rossiter  mode,  or  possibly  higher  harmonics  of  a  lower  Rossiter  mode. 
(It  is  not  possible  to  obtain  accurate  spectra  for  the  higher  frequencies  with  the 
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Figure  5.3:  Mean  from  run  L2.  Each  plot  has  14  contours  equally  spaced  with 
u  G  [-0.2. 0.61],  r  G  [—0.2, 0.2],  h  G  [2.49,2.57],  andct^’  G  [—5,3].  Negative  contours 
are  dashed. 


Mode  # 

9c  energ>'  (u) 

%  energy  (r) 

%  energy  (/i) 

%  energy  (u;) 

1 

33.41 

32.94 

42.76 

31.35 

2 

31.04 

31.09 

29.48 

20.38 

3 

8.62 

9.32 

7.55 

11.92 

4 

8.29 

9.25 

5.49 

9.35 

Table  5.3:  Fraction  of  energ}*  captured  by  each  POD  mode,  run  H2. 
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Figure  5.4:  Mode  1  from  run  L2,  Each  plot  has  14  contours  equally  spaced  with 
u  €  [-1.5. 1.5^.  V  E  [-1.5, 1.5],  h  E  [-1,1],  i*.*  E  -1.5, 1.5],  with  negative  contours 
dashed. 
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Figure  5.6:  Mode  3  from  run  L2.  Contours  as  in  figure  5.4. 
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Figure  5.8:  Eigenvalues  from  POD  of  run  H2:  u  (o),  v  {□),  h  (0),  and  u> 
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Figure  5.9:  Mean  from  run  H2.  Each  plot  has  14  contours  equally  spaced  with 
u  e  [—0.2.0.61],  V  €  [—0.2, 0.2],  h  e  [2.47,2.75],  andcj  G  [—5,3].  Negative  contours 
are  dashed. 

coarse  sampling  used  for  the  snapshots,  listed  in  table  5.1.)  In  general,  higher 
POD  modes  contain  smaller  and  smaller  spatial  scales.  As  before,  we  notice  the 
same  flow  structures  appearing  in  pairs  of  POD  modes,  shifted  in  phase  by  90°. 

Scaling  with  cavity  length 

It  is  interesting  to  compare  (qualitatively)  POD  modes  for  different  operating 
conditions,  to  evaluate  the  scaling  ideas  mentioned  in  section  4.4.4.  For  instance, 
do  the  flow  structures  scale  mainly  with  L/D  or  L/Oq?  We  would  expect  shear 
layer  structures  to  scale  most  strongly  with  L/Oq,  since  shear  layer  characteristics 
do  not  depend  strongly  on  the  cavity  depth  if  the  boundary  layer  is  thin  enough. 

Figure  5.14  supports  this  idea,  that  L/D  is  relatively  unimportant,  compared 
to  L/Oq.  Shown  is  the  first  u-mode  from  run  L2,  next  to  the  corresponding  mode 
from  run  TK4b,  which  is  has  nearly  the  same  value  of  L/6q  and  Reo^  but  with 
L/D  =  4  instead  of  2,  The  pictures  are  indeed  qualitatively  similar,  especially  in 
the  shear-layer  region,  where  the  dynamics  are  important. 

This  scaling  suggests  that  an  appropriate  scaling  for  the  variables  in  the  gov¬ 
erning  equations  described  in  section  4.4.4  is  to  scale  x  by  L  and  y  hy  6o,  thus 
introducing  the  parameter  L/Oq  into  the  equations.  When  nondimensionalized  in 
this  way,  POD  modes  from  runs  with  different  L/D  will  approximately  “collapse 
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onto  each  other,”  for  instance  if  one  simply  ignores  the  bottom  portion  of  the 
L/D  —  2  cavity  in  figure  5.14.  Changing  L/D  will  of  comrse  significantly  change 
the  resonant  acoustic  characteristics  of  the  cavity,  but  the  structures  in  the  shear 
layer  are  more  important  for  capturing  the  correct  d5Tiamics. 

Linear  stability  modes 

The  POD  modes  also  look  qualitatively  similar  to  the  linear  stability  eigenfunctions 
discussed  in  section  3.2.4,  although  these  eigenfunctions  are  obviously  not  orthog¬ 
onal,  Figure  5.15  shows  a  comparison  three  different  types  of  modes,  all  of  which 
look  qualitatively  similar:  POD  modes  1  and  3  from  run  L2;  modes  computed 
from  the  locally-parallel  linear  stability  calculation,  for  the  two  primary  frequen¬ 
cies  observed,  those  of  Rossiter  modes  1  and  2  (St  =  0.4  and  0.7,  respectively); 
and  a  discrete  Fourier  transform  of  the  DNS  data,  for  the  same  two  frequen^-ies' 
(Here  we  have  used  a  tanh  velocity  profile  for  the  linear  stability  calculation.) 

5.1.2  Galerkin  models 

In  this  section,  we  discuss  low-order  models  for  the  cavity,  with  POD  modes  taken 
from  run  H2.  The  equations  we  solve  are  given  in  the  previous  chapter  by  (4.36). 
The  initial  conditions  for  all  runs  are  obtained  by  projecting  a  snapshot  from  the 
D.NS  onto  the  POD  modes,  and  in  this  section  we  always  use  the  snapshot  at 
timestep  70,000,  the  first  in  the  ensemble  used  to  determine  the  POD  modes. 

Effects  of  higher  modes 

Figure  5.16  shows  the  time  coefficients  of  the  first  four  POD  modes,  for  the  projec¬ 
tion  of  the  DNS.  and  for  several  Galerkin  simulations,  where  the  number  of  modes 
varies  from  2  to  20.  Note  that  these  correspond  to  ODEs  with  6  to  60  states,  since 
modes  for  u,  v,  and  h  are  taken  separateljr  Note  that,  as  mentioned  earlier,  the 
time  traces  are  approximately  sinusoidal,  with  the  frequency  of  the  first  two  modes 
corresponding  to  the  first  Rossiter  frequency,  and  the  frequency  of  modes  3  and  4 
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Figure  5.15:  Vorticity  eigenfunctions  computed  from  linear  stability  theory,  for 
first  2  Rossiter  modes  (St  =  0.4. 0.7);  discrete  Fourier  transform  of  DNS  for  the 
same  frequencies;  and  comparison  with  vorticity  POD  modes  1  and  3  for  run  L2. 
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corresponding  to  the  second  Rossiter  frequency.  All  of  the  Galerkin  models  track 
the  DNS  well  for  the  first  three  or  four  periods.  As  we  expect,  for  short  times  the 
20-mode  case  follows  the  DNS  the  closest  for  the  first  few  periods,  but  after  6  or 
7  periods,  all  of  the  models  begin  to  deviate. 

The  same  simulations  are  shown  for  longer  time  in  figure  5.17,  which  demon¬ 
strates  disastrous  long-time  results.  The  4-mode  model  grows  to  an  amplitude 
which  is  much  too  large,  the  10-mode  model  eventually  blows  up,  and  the  20- 
mode  model  demonstrates  complicated,  chaotic-looking  behavior  which  bears  no 
resemblence  to  the  full  DNS  solution.  Surprisingly,  the  2-mode  model  performs 
best  of  all  for  long  time,  stiU  growing  to  an  amplitude  which  is  too  large,  but  at 
least  remaining  stable,  with  the  same  quahtative  features  as  the  DNS  solution. 

The  alarming  result  is  that  increasing  the  number  of  modes  does  not  necessarily 
improve  the  long-time  accuracy  of  the  simulations,  and  can  even  lead  to  instability. 
A  well-known  limitation  of  POD/Galerkin  systems  is  that  while  local  dynamics  are 
guaranteed  to  improve  by  increasing  the  number  of  modes,  there  are  no  guarantees 
about  the  global  dynamical  behavior,  or  even  the  stability  of  equilibrium  points 
and  limit  cycles,  even  if  the  POD  subspace  completely  contains  the  space  spanned 
by  the  snapshots  [Rempfer,  2000]. 

Effects  of  viscosity 

Since  many  of  the  important  phenomena  of  cavity  oscillations  are  essentially  in- 
viscid  (e.g..  shear  layer  instability,  acoustic  wave  propagation),  one  might  expect 
that  the  inviscid  equations  {v  =  0)  would  model  the  flow  nearly  as  well  as  the 
viscous  equations. 

Figure  5.18  shows  time  traces  of  POD  modes  1  and  3,  for  two  different  4-mode 
Galerkin  models,  with  and  without  the  viscous  terms  in  equation  (4.36).  The 
viscous  model  uses  the  same  value  of  Reg  used  in  the  original  DNS  run. 

Both  models  follow  the  DNS  closely  for  the  first  two  periods  of  oscillation,  but 
then  the  amplitude  starts  to  increase,  and  the  inviscid  model  starts  to  perform 
worse.  For  longer  times,  the  inviscid  calculation  eventually  blows  up,  while  the 
viscous  model  remains  stable,  at  least  until  t  =  130,  with  no  sign  of  blow-up. 
For  the  remainder  of  this  thesis,  we  consider  viscous  models  exclusive!}',  with  the 
viscosity  i/  determined  from  Reg  of  the  original  run. 


5.2  Vector- valued  modes 

In  this  section,  we  apply  the  methods  described  in  section  4.4.4  for  obtaining 
reduced-order  models  using  vector-valued  POD  modes. 


5.2  Vector-valued  inodes 


Time  tU/L 

Figure  5.16:  Coefficient  of  u-velocity  modes  1-4,  for  projection  of  DNS  (o),  and 

POD/Galerkin  models  with  2  modes  ( . ),  4  modes  ( - ),  10  modes  ( - ), 

and  20  modes  ( - ). 
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Time  tU/L 


Figure  5.17:  Coefficient  of  u-velocity  modes  1-4,  same  run  as  in  figure  5.1G,  for 
longer  time.  (See  figure  5.16  for  legend.) 
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Figure  5.18:  Coefficient  of  u-velocity  inodes  1  &  3:  4-mode  Galerkin  simulations 
with  viscosity  ( - ),  without  viscosity  ( - ),  and  projection  of  DNS  (o). 

5.2.1  POD  modes 

We  compute  vector- valued  POD  modes  for  the  vector  of  flow  variables  q  =  (u,  v.  a), 
obtaining  the  expansion 


n 

q(x,  y,  t)  =  q(x,  y)  +  qk{t)qkix,  y),  (5.3) 

*=1 

where  q^.  are  the  POD  modes,  and  now  there  is  only  one  set  of  time  coefficients 
(In  the  previous  section^  there  was  a  separate  set  of  time  coefficients  for  each  flow 
variable.)  As  before,  q(x,?/)  is  taken  to  be  the  mean  of  the  snapshots  used  for  the 
POD,  though  this  choice  is  arbitrary. 

In  this  section,  we  do  not  interpolate  onto  a  coarse  grid  first,  but  use  the 
original  grid  from  the  DNS.  However,  we  use  a  truncated  domain  for  computing 
inner  products,  in  particular  the  domain  Q  =  [0,2]  x  [-1,0]  U  [-1,3]  x  [0,2],  the 
same  region  used  in  the  previous  section.  We  use  the  family  of  inner  products 
described  in  section  4.4.3,  with  q  =  1  unless  otherwise  stated  (hence  using  the 
integral  of  the  stagnation  enthalpy  as  a  norm).  Because  the  stretched  grid  is  used, 
the  grid  metrics  must  be  included  in  the  sum,  so  the  inner  products  are  computed 
as 


(qi.q?)  =  ^  {uiU2  +  viV2  + 


2a  ,  dx 


idr} 


A^At], 


(5.4) 
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Figure  5.19:  Normal  velocity  at  a  point  in  the  shear  layer,  for  run  H2;  (o)  indicates 
where  snapshots  are  taken  for  transient  case,  (a)  for  fully-developed  case. 

where  are  coordinates  of  the  uniform  computational  grid,  and  dx/d^,  dy/d-q 
arc  the  grid  metrics.  (Here,  subscripts  i  and  j  denote  evaluation  at  gridpoint 
locations  x,  and  y^.) 

We  consider  two  different  ensembles  of  snapshots  for  computing  the  POD 
modes,  shown  in  figure  5.19,  both  taken  from  run  H2.  For  the  first  set,  which 
we  refer  to  as  the  “transient’’  case,  we  take  snapshots  while  the  flow  is  develop¬ 
ing.  using  201  snapshots  between  0  and  40,000  timesteps.  The  second  set  will  be 
referred  to  as  “fully-developed.”  and  we  take  51  snapshots  between  70.000  and 
80.000  timesteps.  the  same  as  used  for  run  H2  in  the  previous  section.  . 

Figure  5.20  shows  the  eigenvalues  of  the  POD  modes  for  the  two  cases.  As 
we  expect,  the  energj'  decay  is  much  slower  for  the  transient  POD  modes,  as 
more  modes  are  necessary  to  capture  the  more  complicated  flow  structin-es  in  the 
developing  region.  As  before,  the  eigenvalues  occur  in  pairs,  but  for  the  transient 
snapshots,  there  are  additional  POD  modes  corresponding  to  flow  structures  which 
arc  not  oscillating  at  a  particular  frequency,  as  can  be  seen  by  looking  at  the  POD 
modes  thcm.selves. 

Figure  5.21  shows  the  first  sbc  POD  modes  for  the  fully-developed  case,  and 
figure  5.22  shows  the  first  six  modes  for  the  transient  case.  Because  the  modes  are 
vector  valued,  we  can  compute  any  flow  quantity  from  them,  and  here  we  plot  the 
vorticity  and  dilatation,  which  approximately  separates  convecting  disturbances 
from  acoustic  waves.  The  fraction  of  fluctuating  energy  (i|q  -  q]|2)  for  each  mode 
is  also  shown. 

For  the  fully-developed  case,  the  POD  modes  clearly  occur  in  pairs,  shifted  in 
phase  by  90°.  as  in  the  previous  section,  and  the  first  six  modes  together  capture 
99. 59^^  of  the  energy  in  the  fluctuations.  For  the  transient  case,  the  first  six 
modes  capture  much  less  of  the  energy  (93.85%),  and  not  all  of  the  modes  appear 
in  pairs.  In  particular,  modes  1  and  4  both  represent  structures  which  arise  as 
the  flow  changes  from  a  boundary  layer  spanning  a  quiescent  cavity  (the  initial 
condition),  into  an  oscillating  cavity  flow,  with  a  steady  vortex.  For  the  transient 
case,  modes  2  and  3  form  a  pair,  as  do  modes  5  and  6.  These  pairs  may  be 
identified  objectively  by  looking  for  modes  with  similar  energ\'  content.  As  in 
the  scalar-valued  case,  higher  modes  contain  smaller  and  smaller  spatial  scales,  at 
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Figure  5.20:  Eigenvalues  of  POD  modes,  for  snapshots  taken  in  transient  regime 
(a),  and  after  flow  has  fully  developed  (o). 


higher  temporal  frequencies,  with  correspondingly  smaller  acoustic  wavelengths. 

The  fraction  of  energy  captured  by  the  POD  modes  does  not  depend  strongly 
on  the  ensemble  of  snapshots  chosen:  when  a  different  ensemble  of  snapshots  (from 
the  same  run)  is  projected  onto  the  POD  modes,  the  percent  energy  captured  is 
almost  the  same,  as  we  expect,  since  the  dominant  flow  structures  remain  the  same 
once  the  initial  transient  has  decayed. 

The  value  of  a  in  (5.4)  was  varied  from  0  to  2.5,  and  the  POD  modes  were 
virtually  identical  for  all  values  of  a.  This  indicates  that  the  energy  in  the  fluc¬ 
tuations  is  dominated  by  kinetic  energy  in  the  flow  variables,  rather  than  internal 
energ>'  in  the  thermodynamic  variables.  The  POD  modes  are  not  shown,  but  we 
investigate  Galerkin  models  using  different  values  of  a  in  the  next  section. 

5.2.2  Galerkin  models 

We  now  form  reduced-order  models  using  the  vector- valued  POD  modes,  using  the 
equations  given  in  the  previous  chapter  by  (4.45).  As  with  the  scalar- valued  modes, 
the  initial  condition  for  the  simulations  was  obtained  by  projecting  a  snapshot  from 
the  DXS  onto  the  modes. 

Fully  developed  modes 

We  first  examine  models  obtained  from  the  modes  given  in  figure  5.21.  Figure  5.23 
shows  the  results  of  a  simulation  using  an  initial  condition  from  timestep  70,000, 
the  first  snapshot  of  the  ensemble  used  to  compute  the  POD  modes,  and  the  same 
as  used  for  the  initial  condition  in  figure  5.17.  Note  that  since  the  POD  modes  are 


(e)  Mode  5  (0.45%) 


(f)  Mode  6  (0.40%) 


Figure  5.21:  Fully  developed  modes;  Vorticity  (top)  and  dilatation  (bottom),  and 
percent  energj-  captured,  for  run  H2,  with  vector-valued  POD  modes,  snapshots 
from  fully  developed  flow. 


(e)  Mode  5  (1,70%) 


(f)  Mode  6  (1.68%) 


Figure  5.22:  Transient  modes:  Vorticity  (top)  and  dilatation  (bottom),  and  per¬ 
cent  energy  captured,  for  run  H2,  with  vector-valued  POD  modes,  snapshots  from 
developing  region. 
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Figure  5.23:  Fully-developed  modes:  Coefficients  of  POD  modes  1,  3,  and  5  Jor 

Galerkin  simulations,  starting  at  timestep  70,000,  with  2  modes  ( - ),  ••), 

0  ( - ).  11  ( - ),  and  20  modes  ( - ), 


different,  one  cannot  compare  the  time  coefficients  directly  between  the  figures, 
hut  several  qualitative  differences  are  evident. 

In  figure  5.23.  all  of  the  models  track  well  for  the  first  three  or  four  periods. 
After  this,  the  2-mode  model  starts  to  diverge,  increasing  in  amplitude,  but  all  of 
the  simulations  with  more  modes  have  good  qualitative  agreement,  even  for  long 
times.  The  amplitude  of  the  limit  cycle  is  accurately  captured,  and  increasing  the 
number  of  modes  does  not  lead  to  long-time  instability,  as  it  did  with  the  scalar¬ 
valued  modes  shown  in  figure  5.17.  Also  note  that  the  number  of  ODEs  to  solve 
here  is  fewer  by  a  factor  of  3  compared  to  the  scalar-valued  method. 

Figure  5.24  shows  the  same  models  started  from  a  different  initial  condition, 
the  initial  snapshot  from  the  DNS  at  i  =  0,  which  is  a  Blasius  boundary  layer 
spanning  a  quiescent  cavity.  Here,  the  short-time  agreement  is  not  good,  but  the 
qualitative  features  are  correct:  oscillations  gradually  increasing  in  amplitude,  and 
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Figure  5.24:  Fully-developed  inodes:  Time  coefficients  for  Galerkin  simulations, 
starting  at  time  ^  =  0.  Legend  as  in  figure  5.23,  except  2-mode  simulation  not 
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Figure  5.25:  Fully-developed  modes:  Time  coefficients,  varying  a  in  the  inner 
product.  For  a  =  0,  projection  of  DNS  (□),  and  6-mode  Galerkin  simulation 
{ - );  for  a  =  2.5,  projection  of  DNS  (a)  and  6-mode  Galerkin  ( - ). 

ultimately  saturating.  It  also  appears  the  frequency  is  correct,  at  least  for  mode  1, 
and  it  is  difficult  to  tell  the  frequency  content  of  the  higher  modes  in  the  DNS, 
because  the  sample  rate  of  the  snapshots  is  relatively  low.  The  20-mode  case  starts 
to  deviate  around  time  t  =  25,  but  it  does  not  blow  up  (at  least  for  the  length 
of  our  simulation,  which  was  t  =  120).  All  of  the  other  cases  capture  the  final 
amplitude  reasonably  well,  and  increasing  the  number  of  modes  from  4  to  6  to  11 
seems  to  improve  the  agreement  for  the  higher  modes. 

In  figure  5.25,  we  vary  the  value  of  q  used  in  the  inner  product  (5.4),  from  0, 
neglecting  the  thermodynamic  variables  completely',  to  2.5,  weighting  the  thermo¬ 
dynamic  variables  more  strongly’.  The  modes  used  for  these  simulations  were  still 
the  POD  modes  shown  in  figure  5.21 — as  noted  earlier,  changing  a  had  negligible 
effect  on  the  POD  modes,  but  the  value  of  a  also  affects  the  projection  of  the 
equations,  which  is  what  we  investigate  here. 
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As  the  figure  shows,  the  projection  of  the  DNS  snapshots  is  v^ry  similar  for  the 
two  inner  products,  as  is  the  evolution  of  the  diflFerent  sets  of  Galerkin  equations 
obtained  from  the  two  different  projections.  As  mentioned  earlier,  this  indicates 
that  the  energy  in  the  fluctuations  is  dominated  by  the  kinematic  variables,  not 
the  thermodynamic  variables,  which  is  to  be  expected  for  a  cold,  low-Mach  number 
flow  such  as  this.  Since  the  results  are  so  weakly  dependent  on  q,  for  the  remainder 
of  the  section  we  fix  a  =  1.4. 

Transient  modes 

We  now  discuss  simulations  using  POD  modes  from  the  transient  snapshots,  shown 
in  figure  5,22.  Figure  5.26  shows  simulations  varying  the  number  of  modes  between 
3  and  15,  starting  at  timestep  70,000.  Here,  increasing  the  number  of  modes  to 
15  causes  the  flow  to  diverge  after  a  few  oscillations,  unlike  for  the  fully-developed 
modes  (figure  5.23).  Nevertheless,  all  of  the  other  simulations  in  figure  5.26  track 
reasonably  well  for  long  times,  and  none  of  the  simulations  blow  up,  as.  they  did 
for  the  scalar- valued  modes  (figure  5.17). 

Figure  5.27  shows  the  same  simulations  started  from  time  i  =  0.  Here,  all 
of  the  simulations  are  again  accurate  for  short  times,  but  both  the  10-mode  and 
15-mode  cases  deviate  for  long  times.  The  3-mode  case  overshoots  in  amplitude 
(both  modes  1  and  2),  but  eventually  settles  into  a  more  realistic  flow  (well  beyond 
the  time  shown  in  the  figure),  and  the  6-mode  case  tracks  very  well  for  the  entire 
duration  of  the  simulation. 

For  the  simulations  performed  here,  using  the  POD  modes  from  the  developing 
region  gives  worse  results  than  using  the  fully-developed  POD  modes:  more  modes 
are  required  for  qualitatively  accurate  results,  and  furthermore,  taking  more  modes 
tends  to  make  the  models  unstable.  This  is  surprising,  as  one  might  think  that 
the  developing  POD  modes  might  make  Galerkin  projections  more  stable:  since 
they  contain  more  information  about  the  transient?;  they  might  be  better  able  to 
correct  if  the  solution  started  to  stray  from  its  steady-state  conditions.  It  is  likely 
that  many  more  modes  are  necessary,  however,  given  the  much  slower  energ>^  decay 
shown  in  figure  5.20.  The  highest  number  of  modes  we  tried  in  our  simulations 
with  transient  POD  modes  was  15,  and  though  the  first  15  transient  POD  modes 
capture  98.57%  of  the  energ>^  in  the  transient  dataset,  they  presumably  capture 
much  less  energ>'  for  the  fully  developed  oscillations. 

5.2.3  Parameter  variation 

Here  we  briefly  discuss  the  effect  of  varying  the  Mach  number  M  in  the  Galerkin 
models.  In  particular,  we  investigate  how  the  frequency  of  oscillation  varies  as  the 
Mach  number  is  changed,  by  examining  the  eigenvalues  of  the  linearized  system, 
which  is  just  the  linear  part  of  the  governing  equations  (4.45).  For  the  models 
in  this  section,  we  use  POD  modes  from  the  fully-developed  snapshots,  shown  in 
figure  5.21. 

The  eigenvalues  of  the  linearized  system  are  shown  in  Figure  5.28,  for  Galerkin 
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Figure  5.26:  Transient  modes:  Coefficients  of  POD  modes  1,  2,  4,  and  5  for 

Galerkin  simulations,  starting  at  timestep  70,000,  with  3  modes  6  ( - ), 

10  ( - ),  and  15  modes  ( - ). 


Mode  5  Mode  4  Mode  2  Mode  1 
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Figure  5.27:  Transient  modes:  Time  coefficients  for  Galerkin  simulations,  starting 
at  time  t  =  0.  Legend  as  in  figure  5.26. 
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Figure  5.28:  Eigenvalues  of  4-mode  (a),  6-mode  (o),  11-mode  (□),  and  20-mode  (o) 
linearized  Galerkin  models,  at  the  nominal  value  M  =  0.6;  and  of  the  11-mode 
model  with  0.3  <  M  <  0.6  ( - ).  and  0.6  <  M  <  0.9  ( - ). 

models  with  4,  6,  11,  and  20  POD  modes.  The  4-  and  6-mode  models  have  one 
pair  of  unstable  (right-half-plane)  eigenvalues,  at  the  frequency  of  Rossiter  mode  II. 
The  10-  and  20-mode  models  have  an  additional  pair  of  unstable  eigenvalues,  at 
the  lower  frequency  of  Rossiter  mode  I.  This  indicates  that  Rossiter  mode  I  is  still 
present  in  the  dataset,  bin  at  a  much  smaller  energy  than  Rossiter  mode  II. 

As  the  Mach  number  is  varied,  the  eigenvalues  move,  as  shown  in  figure  5.28. 
When  the  Mach  number  is  decreased,  the  eigenvalues  become  more  stable,  w'hich 
is  consistent  with  observations  from  the  simulations. 

Figure  5.29  shows  how  the  frequencies  (imaginary  parts  of  the  unstable  eigen¬ 
values)  change  as  the  Mach  number  is  varied.  Also  plotted  are  the  frequencies 
measured  from  the  DNS  at  M  —  0.6.  and  the  curves  of  the  first  two  frequencies 
predicted  by  Rossiter’s  formula  (3.2).  using  Rossiter’s  original  values  7  =  0.25, 
1/k  =  1.75.  The  frequencies  closely  match  the  frequency  measured  in  the  DNS, 
and  capture  the  same  trend  with  Mach  number  predicted  by  Rossiter’s  model. 
Note  that  we  do  not  expect  our  Galerkin  model  to  be  valid  very  far  from  the 
parameter  range  where  the  snapshots  were  taken,  since  for  different  parameter 
values,  the  flow  structures  will  presumably  be  different. 


5.3  Conclusions 

In  this  chapter,  we  have  compared  reduced-order  models  obtained  from  POD  and 
Galerkin  projection,  and  in  particular  we  have  compared  the  procedure  based 
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Figure  5.29:  Frequencies  predicted  by  Galerkin  model  as  Mach  number  is  varied: 
Im(A/27r)  for  unstable  eigenvalues  A  of  11-mode  Galerkin  model  ( - );  frequen¬ 

cies  measured  from  DNS  (x):  and  frequencies  predicted  by  Rossiter’s  equation 
(— ). 

on  scalar-valued  POD  modes  to  that  based  on  vector-valued  POD  modes.  Both 
methods  capture  the  d\Tiamics  well  for  short  time,  but  the  scalar-valued  models 
in  particular  tend  to  diverge  after  a  few  oscillations.  In  addition,  the  scalar-valued 
modes  tend  to  perform  worse  when  more  POD  modes  are  retained,  in  that  the 
long-time  behavior  is  either  non-physical  or  leads  to  bWv-up. 

A  possible  reason  why  the  scalar- valued  modes  perform  poorly  for  long  time  is 
the  following.  One  of  the  terms  appearing  in  the  governing  equations  is  the  dilata¬ 
tion  For  low-Mach  number  flows  such  as  this,  the  dilatation  remains  small, 

even  though  the  individual  terms  and  Vy  are  not  small.  If  scalar- valued  modes 
are  used,  u  and  i’  may  drift  apart  slightly  after  some  time,  so  that  they  no  longer 
cancel  each  other,  and  the  resulting  error  in  dilatation  can  drive  the  equations  in  a 
non-physical  way.  When  vector- valued  modes  are  used,  the  dilatation  is  computed 
much  more  accurately,  as  it  is  effectively  computed  for  each  POD  mode,  and  not 
from  a  sum  of  POD  modes  for  different  variables  which  may  drift  apart. 

In  order  to  use  these  models  for  control  analysis  or  synthesis,  it  is  of  course 
necessary  to  introduce  the  effects  of  actuation  into  the  models.  The  presence 
of  an  actuator  will  presumably  change  the  flow  structures,  so  POD  modes  need 
to  be  taiken  in  the  actuated  flow,  perhaps  stacking  snapshots  from  different  runs. 
Precisely  how  to  do  this  remains  an  open  question,  and  is  a  topic  for  future  research, 
but  the  present  results  are  promising:  since  the  unactuated  flow  can  be  accurately 
modeled  by  as  few  as  3  POD  modes,  it  is  likely  the  actuated  flow  may  also  be 
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described  by  a  model  of  very  low  dimension. 


Chapter  6 


Control  Design  and  Experimental  Results 


As  mentioned  in  the  introduction,  previous  experiments  at  controlling  cavity  oscil¬ 
lations  have  met  with  limited  success.  For  instance,  a  control  system  might  reduce 
till'  steadv-state  amplitude  at  one  resonant  frequency,  but  increase  the  amplitude 
at  another  frequency  [Williams  &  Fabris,  2000b].  The  goal  of  this  chapter  is  to 
understand  these  effects  using  physics-based  models,  and  to  use  these  models  both 
to  design  control  laws,  and  to  understand  any  performance  limitations. 

We  present  a  linear  model  of  the  cavity,  designed  around  the  experiments  of 
Williams  1'  Fabris  [2000a|.  We  then  design  an  [Hi)  optimal  controller  based  on 
thi.s  model,  and  implement  several  different  control  laws  on  the  experiment.  The 
model  e.xplains  several  phenomena  observed  in  the  experiments,  as  well  as  some 
p<‘rformance  limitations,  discussed  in  section  6.4.3. 

The  results  described  in  this  chapter  represent  joint  work  with  Dave  Williams 
at  Illinois  Institute  of  Technology,  and  we  mention  that  the  experimental  results 
are  vr^-v  recent,  and  have  not  been  fully  analyzed.  Nevertheless,  the  preliminary 
reMil's  rove.d  interesting  features  of  the  controlled  flow,  as  we  discuss  in  section  6.4. 

6.1  Experimental  setup 

The  exitorimental  apparatus  described  in  this  section  was  designed  and  built  by 
David  Williams,  Drazen  Fabris,  and  others  at  the  Fluid  Dynamics  Research  Center 
at  Illinois  Institute  of  Technology.  Experiments  were  performed  using  the  3  ft  x  3  ft 
subsonic  wind  tunnel  at  the  United  States  Air  Force  Academy  in  Colorado  Springs. 
We  give  a  brief  description  of  the  relevant  aspects  of  the  experiment  here,  and  refer 
to  Williams  k.  Fabris  [2000a]  for  further  details. 

A  cavity  model  20  in  long,  15  in  wide,  and  4  in  deep  was  installed  in  the  floor 
of  the  3  ft  X  3  ft  test  section  of  the  wind  tunnel.  With  the  cavity  in  place,  the 
tunnel  has  a  Mach  number  range  of  0.2-0.55.  A  photograph  of  the  cavity  interior 
is  shown  in  figure  6,1,  and  a  photograph  of  the  assembly  from  outside  the  tunnel 
is  shown  in  figure  6.2. 


105 


106 


6  Control  Design  and  Experimental  Results 


Fi,iz:uro  6.1:  Photograph  of  cavity,  from  downstream  and  above. 


Figure  6.2;  Photograph  of  cavity  assembly  from  outside  wind  tunnel 


6.1  Experimental  setup 
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Figure  6.3:  Diagram  of  experimental  apparatus  (side  view).  Location  of  Kulite 
pressure  transducers  is  indicated  by  K1-K8,  location  of  hot  films  indicated  by 
H\VL  HW2. 

6.1.1  Sensors 

The  cavity  was  instrumented  with  eight  Kuhte  pressure  transducers  placed  along 
the  cavity  walls,  one  on  the  upstream  wall,  one  on  the  downstream  wall,  and 
six  along  the  floor,  approximately  equally  spaced.  The  approximate  locations  are 
indicated  in  figure  6.3. 

Velocity  measurements  were  obtained  using  two  hot-film  sensors  placed  in  the 
shear  layer  spanning  the  cavity,  one  near  the  upstream  corner  (Vs in  from  the  wall), 
and  one  near  the  downstream  corner  (liV32in  from  the  wall),  both  in  line  with 
the  lip  of  the  cavity,  and  close  to  the  center  in  the  spanwise  direction.  These  are 
indicated  in  figure  6.3,  and  also  appear  in  the  photograph  in  figure  6.1. 

All  signals  were  passed  through  anti-aliasing  filters  prior  to  sampling  by  a 
digital  data  acquisition  system.  Data  were  sampled  at  6  kHz,  typically  for  65,536 
samples  (10.9 sec).  The  anti-aliasing  filters  were  4th-order  Butterworth  bandpass 
filters,  with  a  pass  band  of  0.6  Hz-2.2  kHz.  (The  high-pass  filter  was  necessary  to 
remove  the  DC  offset,  and  when  needed,  the  DC  component  was  meaisured  using 
a  digital  multimeter.) 

6.1.2  Actuator 

The  flow  was  forced  using  zero-net-mass  blowing  through  a  slot  in  the  upstream 
wall  of  the  cavity,  shown  in  figure  6.3.  The  actuator  was  a  pair  of  500- Watt  Sin 
diameter  loudspeakers  in  an  enclosed  chamber  (visible  in  figure  6.2).  The  loud¬ 
speakers  were  powered  by  a  300-watt  power  amplifier.  Though  the  actuator  injects 
zero  net  mass  through  the  slot,  a  nonzero  net  momentum  is  induced  by  spanwise 
vortices  generated  by  periodically  blowing  through  the  slot  (the  “synthetic  jet” 
effect). 

Earlier  experiments  by  Williams  Fabris  [2000a]  investigated  the  effects  of 
different  orientations  of  the  slot,  and  found  that  the  actuator  was  more  effective 
when  the  slot  was  oriented  in  the  streamwise  direction  (as  indicated  in  the  figure). 
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than  when  it  was  oriented  in  the  wall-normal  direction,  or  at  a  45°  angle  relative 
to  the  free-stream.  Note,  however,  that  Shaw  [1998]  found  the  opposite  trend  for 
a  pulsed-blowing  actuator  (positive  net  mass  injection):  forcing  in  the  w’all- normal 
direction  was  most  effective,  and  forcing  in  the  streamwise  direction  was  least 
effective. 

Though  it  is  possible  that  the  presence  of  the  actuator  assembly  may  introduce 
other  effects  into  the  flow,  for  instance  by  changing  the  structural  resonance  of  the 
cavity,  these  extraneous  effects  are  small,  as  the  model  has  been  tested  with  and 
without  the  actuator  assembly  in  place,  and  the  observed  results  are  identical. 

6.1.3’  Control  implementation 

Both  analog  and  digital  controllers  were  implemented.  The  analog  controller  con¬ 
sisted  of  a  bandpass  filter  and  a  phase  shifter.  The  gain  and  phase  could  be 
continuously  adjusted,  and  the  frequencies  of  the  passbands  could  be  adjusted  in 
discrete  increments. 

Digital  controllers  were  implemented  using  a  dSPACE  interface  board,  running 
on  a  separate  computer  from  the  data  acquisition  system.  For  typical  controllers  we 
were  running,  the  maximum  sample  rate  of  the  dSPACE  system  was  about  20  kHz. 
The  digital  controllers  implemented  will  be  described  in  sections  6.3  and  6.4. 


6.2  Linear  model 

Our  model  for  the  cavity  dynamics  is  based  on  the  familiar  Rossiter  mechanism 
described  in  the  introduction:  small  disturbances  to  the  shear  layer  excite  Kelvin- 
Helmholtz  instabilities,  and  grow  as  they  convect  downstream.  When  these  distur¬ 
bances  impinge  on  the  trailing  edge  of  the  cavity,  they  are  scattered  into  acoustic 
waves,  which  propagate  upstream  and  excite  further  instabilities  in  the  shear  layer. 
This  process  is  depicted  in  the  block  diagram  shown  in  figure  6.4,  where  we  have 
written  each  component  of  the  mechanism  described  above  as  a  separate  transfer 
function. 

Here.  G{s)  represents  the  shear-layer  transfer  function,  i.e.,  the  transfer  func¬ 
tion  from  \'elocity  disturbances  uq  at  the  leading  edge  to  velocitv  disturbances 
at  the  trailing  edge.  Transfer  functions  for  acoustic  scattering,  propagation,  and 
recepti\  it\  are  given  by  S,  A,  and  R,  and  in  the  diagram  po  3.nd  pi^  denote  pressure 
disturbances  at  the  leading  and  trailing  edges,  respectively.  These  quantities  may 
be  measured  from  the  experiment:  vo  is  measured  by  hot  film  1,  wl  by  hot  film  2, 
and  Pi  and  po  by  Kulites  2  and  8,  respectively  (see  figure  6.3).  Here,  we  do  not 
use  Kulite  1,  as  this  sensor  measures  substantial  pressure  fluctuations  from  the 
impinging  shear  layer. 

The  other  transfer  functions  depicted  in  figure  6.4  represent  the  influence  of 
the  actuator  and  controller.  The  control  output  u  is  the  voltage  to  the  amplifier, 
and  we  use  the  signal  from  Kulite  8  as  the  input  y.  The  controller  transfer  function 
(which  we  choose)  is  given  by  C,  and  the  actuator  dynamics  are  described  by  a 
transfer  function  V. 
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Figure  6.4:  Block  diagram  of  cavity  model. 
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Figure  6.5:  Block  diagram  of  cavity  model  with  disturbances. 


The  overall  transfer  function  P  for  the  cavity  is  then 

(6.1) 

For  the  purposes  of  studying  the  dvmamical  features  of  this  model,  we  ignore  the 
actuator  d\Tiamics,  setting  V  —  1.  (These  actuator  dynamics  may  be  measured 
from  the  experiment,  and  once  measured,  their  effects  may  be  inverted  out  of  the 
control  laws  we  obtain.)  Theoretical  models  for  the  remaining  transfer  functions 
are  discussed  below. 

Figure  6.4  does  not  include  the  effects  of  disturbances.  Our  model  for  the  way 
in  which  disturbances  enter  is  the  standard  approach  shown  in  figure  6.5,  where 
n  is  the  sensor  noise,  and  d  is  the  process  noise.  While  it  might  make  better 
physical  sense  for  the  process  noise  to  enter  after  the  actuator  dynamics  V  (e.g., 
as  a  disturbance  to  vo  in  figure  6.4),  this  effect  may  easily  be  inverted  by  changing 
the  frequency  content  assumed  for  d. 
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x/L  =  0.0  x/L  =  0.1  x/L  =  0.2  x/L  =  0.3  x/L  =  0.4  x/L  =  0  5  x/L  =  0.6 


Figure  6.6:  Velocity  profiles  for  the  cavity  shear  layer.  Hot  wire  measurements  (□) 
and  tanh  profiles  with  same  vorticity  thickness  and  deflection  ( - ). 

6.2.1  Shear  layer 

The  shear  layer  transfer  function  G{s)  is  determined  from  linear  stability  theory. 
We  begin  with  velocity  profiles  measured  in  experiments  by  Williams  Fabris 
[2000a].  shown  in  figure  6.6.  These  profiles  are  from  a  run  with  Mach  number  M  = 
0.35.  in  a  cavity  with  aspect  ratio  L/D  =  5.  Figure  6.6  show's  the  experimental 
data  along  w-ith  hyperbolic  tangent  profiles  with  the  same  vorticity  thickness.  The 
spreading  rate  of  the  shear  layer  is  determined  from  a  linear  fit  to  the  data,  and 
used  as  an  input  to  a  hnear  stability  calculation  to  determine  the  amplification  and 
phase  of  shear  layer  disturbances.  We  then  fit  a  rational  function  to  the  resulting 
transfer  function,  and  the  result  is  shown  in  figure  6.7. 

6.2.2  Acoustics 

The  model  we  use  for  acoustic  propagation  in  the  cavity  is  show'n  in  figure  6.8. 
Here,  r  =  L/a  is  a  time  delay  which  represents  the  acoustic  lag  between  the  trailing 
edge  and  leading  edge  (here,  L  is  the  cavity  length  and  a  is  the  sound  speed  inside 
the  cavity).  An  acoustic  wave  emanating  from  the  downstream  corner  x  =  L 
propagates  upstream,  and  some  of  it  reflects  off  the  upstream  wall,  propagates 
downstream,  and  again  reflects  off  the  downstream  wall.  The  reflection  coefficient 
r  measures  the  attenuation  in  these  reflections  (e.g.,  if  both  reflections  are  perfect, 
then  r  =  1;  if  each  reflection  reduces  the  amplitude  by  0.5,  then  r  =  0.25).  This 
model  therefore  captures  longitudinal  modes  of  acoustic  resonance,  but  ignores 
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Figure  6.7:  Bode  plot  of  shear  layer  transfer  function  G^is),  determined  from  linear 
stability  theory. 
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Figure  6.8:  Block  diagram  of  transfer  function  >l(s)  for  cavity  acoustics. 


112 


6  Control  Design  and  Experimental  Results 


Figure  6.9:  Frequency  response  of  acoustic  transfer  function,  for  M  =  0.6.  Reflec¬ 
tion  coefficient  r  =  0  ( - );  r  =  0.25  ( - );  r  =  0.5  (— );  r  =  0.75  ( . ). 

Note  the  linear  frequency  scaJe. 

depth  modes.  It  is  probably  reasonable  to  ignore  the  depth  modes  for  such  a 
shallow  cavity  {L/D  =  5). 

The  overall  transfer  function  is  given  by 


and  the  Bode  plot  is  shown  in  figure  6.9,  for  M  =  U/a  =  0.6,  and  various  values 
of  r.  ranging  from  0  to  0.75.  The  resonant  peaks  are  clearly  apparent,  and  for 
this  simple  model,  all  of  the  harmonics  are  equally  strong.  Note  that  these  res¬ 
onant  peaks  represent  the  longitudinal  acoustic  modes  of  the  cavity,  and  not  the 
Rossiter  frequencies.  However,  when  these  resonant  acoustic  frequencies  approach 
the  Rossiter  frequencies,  they  may  influence  the  mode  selection,  determining  which 
Ro.ssiter  mode  is  dominant.  ° 

6.2.3  Scattering  and  Receptivity 

The  scattering  and  receptivity  effects,  which  couple  vortical  and  acoustic  distur¬ 
bances  at  the  trailing  and  leading  edge,  are  the  least  well  understood  aspects  of 
the  cavity  model.  Details  of  these  effects  have  been  studied  by  Crighton  [1992] 
for  edge  tones,  and  Kerschen  [1996]  for  cavity  flows.  However,  both  of  these  mod¬ 
els  describe  the  scattering  by  a  sharp  edge,  rather  than  a  corner;  scattering  at  a 
sharp  edge  produces  a  dipole  source,  as  is  well  known  for  edge  tones,  while  the 
acoustic  source  in  the  cavity  is  more  closely  represented  by  a  monopole,  as  men¬ 
tioned  in  section  3.4.3.  Furthermore,  these  previous  scattering  models  are  quite 
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detailed,  employing  a  Wiener-Hopf  factorization  (which  does  not  extend  easily  to 
corners),  and  for  the  purposes  of  control  we  are  not  concerned  with  detailed  flow 
features,  but  merely  the  overall  phase  and  amphtude  effects,  as  a  function  of  fre¬ 
quency.  Hence,  in  this  section  we  consider  a  crude  model,  based  on  Rossiter’s 
semi-empirical  formula.  .  ^ 

Rossiter's  formula,  given  by  equation  (3.2),  is  derived  as  follows.  Consider  the 
total  phase  variation  tp  around  the  cavity  feedback  loop,  given  by 

ip  =  aL  —  27r7  -I-  kL,  (6.2) 

where  a  is  the  wavenumber  of  convected  disturbances  in  the  shear  layer,  k  is  the 
wavenumber  of  acoustic  waves,  and  7  is  an  empirical  constant  which  represents 
a  phase  lag,  which  Rossiter  attributed  to  receptivity  and  scattering  effects.  The 
wavenumbers  may  be  rewritten  in  terms  of  the  shear  layer  convection  speed  Cp  = 
— a;/Q  and  the  sound  speed  a  =  —uj/k.  Resonant  frequencies  occur  whenever 
=  -27rn,  so  solving  (6.2)  for  uj,  we  obtain  Rossiter’s  formula  for  the  resonant 
frequencies 

^  n-7 

2nU  M+1/ac’ 

where  k  —  Cp/U.  Rossiter  found  that  with  7  =  0.25,  a  good  correlation  was  ob¬ 
tained  between  the  resonant  frequencies  given  this  formula  and  those  measured  in 
experiment.  As  shown  by  (6.2),  the  corresponding  phase  shift  from  scattering  and 
receptivity  is  ~-7r/2,  independent  of  frequency.  The  analytic  transfer  function  S{s) 
with  this  phase  shift  is  simply  an  integrator,  S(s)  =  1/s.  This  transfer  function 
has  an  infinite  gain  at  zero  frequency,  which  is  physically  unreasonable,  but  this 
may  be  avoided  by  using  a  first-order  lag  instead  of  a  pure  integrator: 

5(5)  =  '  (6.4) 

s  4-  a 

where  A’  is  a  gain,  and  a  is  cutoff  frequency  which  is  small  compared  to  all  resonant 
frequencies. 

For  the  present  mode,  we  use  (6.4)  for  the  scattering  effects,  and  for  the  recep¬ 
tivity  we  take  B  =  1.  This  is  undoubtedly  a  crude  model,  but  no  more  crude  than 
Rossiter's  model,  which  has  been  shown  to  predict  resonant  frequencies  reasonably 
well. 

6.2.4  Overall  cavity  model 

The  overall  cavity  transfer  function  P  is  then  formed  from  equation  (6.1),  and  its 
poles  and  zeros  are  plotted  in  figure  6.10.  (For  the  time  delay,  we  use  a  6th-order 
Fade  approximation  to  obtain  a  rational  transfer  function.)  The  imaginary  part 
of  the  poles  neetr  the  imaginary  axis  closely  match  the  Rossiter  frequencies,  as  we 
expect.  The  first  three  Rossiter  frequencies  are  unstable,  and  higher  frequencies 
are  stable,  since  they  are  damped  by  the  shear  layer.  The  unstable  frequencies 
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Figure  6.10:  Poles  (x)  and  zeros  (o)  of  the  cavity  transfer  function  P{s).  A  6th- 
order  Fade  approximation  is  used  for  the  time  delay.  Note  unstable  poles  at  the 
first  three  Rossiter  frequencies. 


predicted  by  our  model  (at  M  =  0.35)  are  100  Hz,  207  Hz,  and  312  Hz.  In  the 
experiment,  onl\  the  second  and  third  modes  were  observed,  at  frecjuencies  of 
about  190  Hz  and  337  Hz. 


6.3  Control  design 

In  this  section,  we  apply  standard  tools  of  control  design  to  the  cavity  model 
developed  above.  We  use  Linear  Quadratic  Gaussian  (LQG)  methods  to  design 
a  controller  which  is  optimal,  based  on  prescribed  quadratic  cost  functions  and 
simple  assumptions  about  the  disturbances. 

6.3.1  LQG  design 

We  begin  with  a  state-space  realization  of  the  cavdty  model  P{s),  written  in  the 
form 


X  =  Ax  +  Bu 
y  -  Cx, 


(6.5) 


where  u  and  y  are  the  inputs  and  outputs,  respectively  (see  figure  6.4),  and  x 
is  the  state  variable.  In  the  LQG  design  procedure,  one  first  designs  a  state- 
feedback  controller  u  -  Kx  which  stabilizes  the  system,  assuming  full  knowledge 
of  the  state,  and  then  one  designs  an  observer  to  estimate  the  state  based  on  output 
measurements  y.  The  controller  gains  K  are  found  by  solving  the  Linear  Quadratic 
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Regulator  (LQR)  problem,  which  gives  a  hneax  controller  which  is  optimal,  based 
on  a  prescribed  quadratic  cost  function.  To  estimate  the  state,  one  uses  a  Kalman 
filter,  which  is  a  state  estimator  which  is  optimal  in  terms  of  minimizing  the  mean 
square  estimation  error,  based  on  some  assumptions  about  the  disturbances. 

The  cost  function  we  minimize  for  the  LQR  control  design  is 

CX) 

(Qy^  +  u^)  dt  (6.6) 

where  Q  is  an  arbitrary  weight  which  determines  how  much  to  penahze  the  control 
input,  relative  to  the  output  error. 

The  state  estimator  provides  an  estimate  x  of  the  exact  state  x  (which  cannot 
be  measured  directly).  The  estimate  is  determined  by 

X  =  Ax  ^  Bu  +  L{y  -  y) 
y  =  Cx. 

a  system  which  mimics  (6.5),  but  is  forced  by  the  output  error  y  ~y.  It  is  simple 
to  show  that  the  error  x  —  x  converges  to  zero  as  long  as  the  eigenvalues  oi  A  — 
LC  are  stable.  The  Kalman  filter  specifies  observer  gains  L  which  are  optimal, 
in  the  sense  that  the  error  converges  the  fastest  in  the  presence  of  stochastic 
disturbances  d  and  n  (see  figure  6.5),  assumed  to  be  Gaussian  white  noise,  with 
specified  covariances. 

The  LQR  and  optimal  estimator  problems  are  easily  solved  using  standard 
routines  in  Matlab.  The  controller  and  estimator  are  then  combined  to  give  a 
compensator,  which  may  be  written  as  a  transfer  function  C{s)  of  the  same  order 
as  the  original  plant  P{s).  This  order  of  the  compensator  may  be  reduced  using 
ba^'-mced  truncation,  or  more  crudely  by  simply  canceling  nearb}*  pole-zero  pairs. 

6.3.2  Stabilizing  controller  for  the  cavity  model 

In  the  model  shown  in  figure  6.4,  it  is  clear  that  a  controller  which  merely  cancels 
the  natural  acoustic  feedback  in  the  cavity  will  stabilize  the  plant.  Any  process 
noise  d  will  still  be  amplified  by  the  shear  layer,  but  the  feedback  path  which  leads 
to  cavity  resonance  will  be  broken. 

If  we  neglect  the  actuator  d>Tiamics  (V'  =  1),  then  a  stabilizing  controller  is 
given  by  C  =  — i?“L  Since  we  take  /?  =  1  in  the  model,  the  stabilizing  controller 
is  just  C  —  —1.  Indeed,  the  optimal  control  techniques  outlined  in  the  previous 
section  give  a  compensator  which  closely  matches  this  result. 

The  optimal  controller  given  by  LQG  is  very  high-order  (the  same  dimension 
as  the  state,  which  is  high-order,  as  can  be  seen  by  the  large  number  of  poles  and 
zeros  in  figure  6.10),  so  we  perform  an  elementary  model  reduction  on  the  compen¬ 
sator  C(s),  removing  (stable)  pole-zero  pairs  which  nearly  cancel  each  other.  After 
model  reduction,  the  compensator  is  sixth-order  (reduced  from  16),  and  its  trans¬ 
fer  function  is  shown  in  figure  6,11.  Note  that,  in  the  frequency  range  where  the 
model  is  valid  (less  than  about  3000rad/sec),  the  compensator  closely  resembles 
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From:  y1  (Measurement) 


Figure  6.11:  Bode  plot  of  controller  transfer  function  C{s)  (reduced-order  con¬ 
troller).  For  reference,  the  unstable  frequencies  occurred  at  uj  —  628,  1301,  and 
1960  rad/s. 


the  function  —1  (180°  phase,  with  a  gain  of  OdB),  as  expected. 

The  poles  of  the  closed-loop  system  (with  control)  are  shown  in  figure  6.12. 
\\  ithout  the  model  reduction  step,  LQG  guarantees  the  closed-loop  system  will 
be  stable,  but  figure  6.12  demonstrates  that  the  closed  loop  system  is  stable  (all 
poles  in  the  left  half-plane)  even  with  the  reduced-order  compensator. 

For  such  a  simple  model,  where  a  stabilizing  control  law  is  obvious,  wc  really 
do  not  require  the  elaborate  machinerj-  of  LQG  design.  Our  methodologv  is  nev¬ 
ertheless  useful,  as  enhancements  to  the  model  (for  instance,  measmiug  transfer 
functions  directly  from  an  experiment)  will  lead  to  more  complex  controllers.  Fur¬ 
thermore,  once  we  have  a  model  of  the  physics,  we  may  analyze  limitations  of 
achierable  performance,  as  we  discuss  in  the  next  section. 


6.4  Experimental  results 

Here,  we  present  some  recent  results  from  the  experiment  described  in  section  6.1, 
discussing  results  using  both  analog  and  digital  controllers.  All  of  the  results  are 
for  a  run  with  M  =  0.34,  at  which  only  a  single  mode  of  cavity  oscillations  was 
present. 

6.4.1  Analog  controller 

The  analog  controller  consisted  of  an  Ithaco  filter,  a  4th-order  Butterworth  band¬ 
pass  filter  whose  passbands  could  be  adjusted  in  discrete  increments,  and  a  phase 
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Figure  6.12:  Poles  and  zeros  of  closed-loop  transfer  function,  with  reduced-order 
controller.  Rossiter  modes  are  stabilized. 


Name 

Order 

Passband 

Filter  1 

2 

320-360  Hz 

Filter  2 

2 

290-390  Hz 

Filter  3 

2 

215-465  Hz 

Table  6.1:  Parameters  of  digital  Butterworth  filters  used. 

shifter,  which  could  be  fine-:  aned  manually.  The  transfer  function  of  the  analog 
controller  is  shown  in  figure  5  13,  where  the  passband  was  set  to  300-350  Hz. 

Figure  6.14  shows  the  results  of  the  experiment,  with  and  without  control. 
For  the  case  showm,  the  gain  and  phase  were  manually  tuned  for  the  greatest 
suppression.  When  control  is  switched  on,  the  resonant  frequency  at  337  Hz  is 
attenuated  significantly,  but  a  new  peak  appears  at  434 Hz.  The  smaller  peah  at 
674  Hz  is  the  first  harmonic  of  the  fundamental  cavity  frequency,  and  it  too  is 
attenuated  with  control. 

6.4.2  Digital  controller 

The  digital  controllers  implemented  were  also  bandpass  filters,  and  we  tried  several 
digital  Butteru^orth  filters  of  different  orders  and  passbands.  The  orders  and  pass- 
bands  for  the  filters  used  here  are  shown  in  table  6.1,  and  the  transfer  functions 
arc  shown  in  figure  6.15. 

Figure  6.16  shows  the  results  of  the  closed-loop  experiments  with  the  different 
filters.  In  this  figure,  the  gain  and  phase  was  tuned  for  the  best  suppression.  The 
narrow  band  filter  showed  very  little  attenuation,  and  the  broadband  filters  better 
attenuated  the  main  cavity  frequency,  but  a  higher  frequency  peah  appears,  as  it 
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Figure  6.13;  Bode  plot  of  analog  controller.  Frequency  of  cavity  oscillations  indi 
cated  by  (□). 


(a)  Signal  from  Kulite  3 


Figure  6.14:  Power  spectrum  from  two  pressure  sensors,  for  analog  controller 
control  off  ( - );  control  on  ( - ). 


6.4  Experimental  results 


Figure  6.15:  Bode  plot  of  filters  used  in  digital  controllers:  filter  1  (* 
( - );  filter  3  ( - ). 
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(a)  Filter  1  (b)  Filter  2  (c)  Filter  3 


Figure '6. 16:  Power  spectrum  with  3  different  digital  filters,  tuning  gain  and  phase 
manually  for  best  suppression:  control  off  ( - );  control  on  ( - ). 


(a)  Delay  =  4  x  10  sec  (b)  Delay  =  5  x  10  sec  (c)  Delay  =  6  x  lO”"*  sec 


Figure  6.17:  Power  spectra  with  digital  controllers,  showing  sidebands:  control  off 
( - ):  control  on,  using  filter  2  ( - ). 

did  with  the  analog  controllers.  The  frequency  of  this  peak  shifts  with  different 
controllers:  with  filter  2,  the  peak  is  at  430 Hz;  with  filter  3,  449 Hz. 

Note  that,  especially  evident  in  figure  6.16b,  the  main  cavity  frequency  is  some¬ 
times  split  into  two  sidebsinds  when  the  control  is  turned  on.  This  effect  is  explored 
further  in  figure  6.17,  where  filter  2  was  used,  and  the  time  delay  was  adjusted. 
The  main  resonant  frequency  at  337  Hz  is  almost  completely  attenuated,  but  side¬ 
bands  appear  very  close  in  frequency,  at  about  320 Hz  and  341  Hz.  As  the  time 
delay  is  changed,  the  relative  strength  of  the  sidebands  changes,  and  the  frequency 
changes  slightly — the  lower  frequency  shifts  from  320  to  325  Hz  in  figures  (a)-(c). 
We  discuss  these  sidebands  further  in  the  next  section. 

6.4.3  Performance  limitations 

The  experimental  results  of  the  previous  section  are  somewhat  discouraging:  when¬ 
ever  significant  attenuation  of  the  main  cavity  frequency  was  achieved,  new  oscilla¬ 
tions  appeared  at  different  frequencies.  These  new  peaks  of  the  closed-loop  system 
are  of  two  t>*pes:  a  broad  peak  occurred  at  a  frequency  quite  different  from  the 
main  cavity  frequency  (e.g.,  at  about  430 Hz  in  figure  6.17);  in  addition,  narrow 
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peaks  very  close  to  the  original  cavity  frequency  occurred,  which  we  refer  to  as 
sidebands. 

A  similar  appearance  of  sidebands  has  been  observed  in  recent  experiments 
on  control  of  combustion  instabilities,  discussed  by  Banaszuk  et  al  [2001].  These 
combustion  instabilities  are  similar  to  cavity  instabilities  in  many  ways:  both  sys¬ 
tems  exhibit  self-sustained  oscillations  caused  by  flow- acoustic  interactions.  Time 
delays  are  inherent  in  both  systems,  because  of  the  nature  of  the  acoustic  teed-  ^ 
back.  In  both  systems,  the  amphtude  of  oscillations  is  determined  by  a  nonlinear 
saturation  mechanism  (though  this  mechanism  is  much  better  understood  in  the 
combustion  problem). 

The  appearance  of  these  sidebands  may  be  understood  from  a  purely  linear 
point  of  view,  employing  the  well-known  “area  rule,”  which  estabhshes  performance 
limitations  for  feedback  systems.  The  area  rule  determines  constraints  on  the 
sensitivity  function 
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1 
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for  linear  systems  of  the  form  6.5,  The  area  rule  applies  to  systems  whose  relative 
degree  is  at  least  2  (relative  degree  is  the  degree  of  the  denominator  of  a  transfer 
function  minus  the  degree  of  the  numerator),  and  states  that 


roc 

logio|S(za;)i(L;  =  Trlogioe^Repi, 


(6.8) 


where  p,  are  the  unstable  poles  of  the  transfer  function  PC.  Thus,  if  disturbances 
are  attenuated  over  some  frequency  range,  they  must  be  amplified  over  some  other 
frequency  range,  and  if  the  plant  is  unstable  (so  that  there  are  unstable  poles  r»i) 
then  this  effect  is  exacerbated.  Furthermore,  this  effect  is  psirticularly  apparent 
for  narrow  bandwidth  controllers.  As  PC  0,  S  1,  so  logjo  |5|  ->  0  Thus, 
all  of  the  disturbance  amplification  must  occur  within  the  narrow  bandwidth  of 
the  controller,  and  the  more  narrow  the  bandwidth,  the  greater  the  amount  of  the 
disturbance  amplification. 

It  Ls  likely  that  this  effect  explains  the  sidebands  evident  in  figure  6.17.  Though 
the  models  developed  in  section  6.2  do  not  represent  the  experiment  closely  enough 
to  verify  this  effect  quantitatively  (our  models  even  predict  more  than  one  fre¬ 
quency  of  oscillation,  while  the  experiment  shows  only  one  mode),  more  refined 
models  may  well  be  able  to  explain  these  phenomena.  It  is  likely  that  more  ac¬ 
curate  models  may  be  obtained  by  using  experiments  to  identify  the  individual 
transfer  functions  represented  in  figure  6.4,  but  we  leave  this  for  future  work. 
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6  Control  Design  and  Experimental  Results 


Chapter  7 
Conclusions 


This  chapter  contains  a  summary  of  the  contributions  of  this  thesis,  along  with 
some  suggestions  for  future  study. 


7.1  Summary 

We  have  performed  an  extensive  set  of  Direct  Numerical  Simulations  (DNS)  for  the 
two-dimensional,  compressible  flow  past  a  rectangular  cavity.  These  simulations 
reveal  two  distinct  modes  of  self-sustained  oscillation:  a  shear-layer  mode,  which 
is  the  mode  originally  described  by  Rossiter  [1964],  and  a  wake  mode,  which  was 
first  classified  by  Gharib  Roshko  [1987].  In  wake  mode,  the  frequency  of  oscil¬ 
lation  is  independent  of  Mach  number,  indicating  that  acoustics  no  longer  play  an 
important  role  in  the  underlying  feedback  mechanism  that  leads  to  self-sustained 
oscillations.  Transition  to  wake  mode  appears  to  occur  when  the  amplitude  of  os¬ 
cillations  exceeds  a  certain  threshold.  In  wake  mode,  there  is  significant  backflow 
in  the  cavity,  and  it  is  possible  that  this  creates  a  region  of  absolute  instabil¬ 
ity  in  the  shear  layer,  which  may  provide  the  feedback  mechanism  necessary  for 
self-sustained  oscillations. 

The  technique  of  Proper  Orthogonal  Decomposition  (POD)  and  Galerkin  pro¬ 
jection  is  well  known,  and  has  been  used  to  study  a  number  of  incompressible  flows 
[Aubry  et  al.,  1988].  We  have  extended  this  technique  to  compressible  flows.  Using 
a  simplified  form  of  the  Navier-Stokes  equations,  valid  for  cold  flows  at  moderate 
Mach  number,  we  have  obtained  Galerkin  projections  which  are  quadratic,  and 
are  easily  computed.  We  have  introduced  a  family  of  inner  products  which  may  be 
used  for  computing  vector-valued  POD  modes  for  compressible  flows,  where  both 
the  kinematic  and  thermodynamic  variables  contribute  to  the  total  energy. 

We  have  applied  these  model-reduction  techniques  to  the  flow  past  a  cavity, 
using  data  from  our  numerical  simulations  to  obtain  POD  modes.  We  compare 
models  obtained  from  separate  sets  of  scalar- valued  POD  modes,  and  a  single  set  of 
vector- valued  POD  modes,  obtaining  reduced-order  systems  with  between  2  and  60 
states.  All  of  the  models  capture  the  dynamics  well  for  short  times  (a  few  periods 
of  the  oscillation),  and  in  general  the  models  based  on  scalar- valued  modes  diverge 
for  long  times,  and  are  even  unstable,  while  the  models  using  vector- valued  modes 
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7  Conclusions 


have  much  better  long-time  behavior. 

Finally,  we  have  described  a  conceptual  modeling  approach,  where  we  model 
the  individual  components  of  the  cavity  physics  separately.  This  approach  leads 
to  linear  models  to  which  standard  tools  from  control  theory  apply,  and  we  use  the 
Linear  Quadratic  Gaussian  (LQG)  technique  to  design  a  controller  which  stabilizes 
cavity  oscillations. 

7.2  Suggestions  for  future  work 

For  the  POD/Galerkin  models  developed  in  this  thesis  to  be  useful  for  bifurcation 
analysis,  the  models  must  first  be  vahdated  over  the  parameter  range  considered. 
To  obtain  low-order  models  valid  over  a  range  of  parameters,  one  might  use  POD 
modes  obtained  from  an  ensemble  of  snapshots  taken  from  several  different  runs 
and  validate  the  resulting  models  by  comparing  to  direct  numerical  simulations 
at  different  parameter  values.  The  resulting  low-order  models  could  then  be  used 
to  study  bifurcations — for  instance,  the  onset  of  cavity  oscillations  qualitatively 
resembles  a  supercritical  Hopf  bifurcation,  and  the  POD/Galerkin  models  could 
be  used  to  provide  a  more  quantitative  understanding  of  this  transition  process. 

Similarly,  for  these  models  to  be  useful  for  control  analysis  or  synthesis,  they 
must  include  the  effects  of  actuation.  One  may  introduce  actuator  effects  using,  for 
instance,  body  force  terms  in  the  equations  of  motion,  but  the  resulting  low-order 
models  also  need  to  be  validated  in  the  presence  of  actuation. 

Finally,  the  experimental  results  presented  in  the  previous  chapter  are  very 
recent,  and  need  to  be  studied  in  greater  detail.  In  particular,  the  experimen¬ 
tal  results  may  be  used  to  identify  the  different  components  of  the  linear  models 
described  in  the  last  chapter,  obtaining  transfer  functions  directly  from  experi¬ 
mental  data.  The  resulting  models  should  rep’^esent  the  experiment  much  more 
faithfully  than  the  simplified  analytical  models  presented,  and  could  be  used  to 
explain  the  effects  observed,  leading  to  a  better  understanding  of  the  flow  physics, 
and  ultimately  guiding  the  design  of  more  effective  control  techniques. 
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